text:
#ifndef FEODOROV_SQC_LAB1_GAUSS_H
#define FEODOROV_SQC_LAB1_GAUSS_H

#include <gsl/gsl_matrix.h>

#define EPSILON .0000001

#ifdef __cplusplus
extern "C" {
#endif

void gauss_init(size_t size, const double *d_matrix, const double *d_vector, gsl_matrix **matrix, gsl_vector **vector, gsl_vector **result);
int gauss(size_t size, gsl_matrix *matrix, gsl_vector *vector, gsl_vector *result, double *det);
void gauss_print(size_t size, gsl_vector *result);
int gauss_check(size_t size, const double *d_matrix, const double *d_vector, gsl_vector *result);

#ifdef __cplusplus
}
#endif

#endif //FEODOROV_SQC_LAB1_GAUSS_H


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

#ifndef rasm_150d55c1_d63d_40ea_8fd2_9bc280975fe0_h
#define rasm_150d55c1_d63d_40ea_8fd2_9bc280975fe0_h

#include <gslib\pool.h>
#include <gslib\string.h>
#include <gslib\std.h>
#include <rathen\config.h>

/*
 * @ 2013.02.26
 * New assemble tool for rathen. This version we gonna solve the problem that last version of rasm could not support labeled jump.
 * In this version, every chuck of assembly code would be a separate and independent node, linked to a whole code segment, and the
 * label would be a particular node in the tree. In this version, we would not need vslen to do the linkage stuff any more.
 * And the short jump would be available.
 * Plus, in this way, we could see the machine code correspond to every single assemble instruction, and yet we've been very close
 * to a disassemble toolkit, but there were still some problems. The rathen didn't use all the x86 machine code. This assemble tool
 * just concerned about a subset of all the x86 machine codes. I didn't intend to provide such a tool to decode all the machine codes
 * into assemble, or rathen source code, either. As the optimize work we might probably did in the future, this decode works would
 * be a hard problem, which was considered unnecessary.
 */

__rathen_begin__

enum reg
{
    _regzr = 0,
    _eax = 1,
    _ecx = 2,
    _edx,
    _ebx,
    _esp,
    _ebp,
    _esi,
    _edi,
    /* more to come... */
};

#ifdef _GS_X86

inline byte get_register_index(reg r)
{
    switch(r)
    {
    case _eax:      return 0;
    case _ecx:      return 1;
    case _edx:      return 2;
    case _ebx:      return 3;
    case _esp:      return 4;   /* why "no index" */
    case _ebp:      return 5;
    case _esi:      return 6;
    case _edi:      return 7;
    default:        assert(0);
    }
    return -1;
}

typedef const gchar*    lab;

#define __rasm_make_type(tname, dname) \
    struct tname \
    { \
        dname       _data; \
        \
        explicit tname(dname d = (dname)0) { _data = d; } \
        bool operator == (const tname& d) const { return _data == d._data; } \
        dname inner_data() const { return _data; } \
    }

__rasm_make_type(_bias_, int);
__rasm_make_type(_bias8_, int8);
__rasm_make_type(_bias16_, int16);
__rasm_make_type(_mem_, int);
__rasm_make_type(_reg_, reg);

enum cccc
{
    cccc_o = 0,
    cccc_no,
    cccc_c,
    cccc_nc,
    cccc_e,
    cccc_ne,
    cccc_be,
    cccc_a,
    cccc_s,
    cccc_ns,
    cccc_p,
    cccc_np,
    cccc_l,
    cccc_ge,
    cccc_le,
    cccc_g,
};

struct __gs_novtable rasm_instr abstract
{
    virtual ~rasm_instr() {}
    virtual const gchar* get_name() const = 0;
    virtual bool finalize(vessel& vsl) = 0;
    virtual void to_string(string& str) const = 0;
    virtual bool compare(const rasm_instr* that) const = 0;
};

inline const gchar* nameofreg(reg r)
{
    switch(r)
    {
    case _eax:  return _t("eax");
    case _ecx:  return _t("ecx");
    case _edx:  return _t("edx");
    case _ebx:  return _t("ebx");
    case _esp:  return _t("esp");
    case _ebp:  return _t("ebp");
    case _esi:  return _t("esi");
    case _edi:  return _t("edi");
    }
    return 0;
}

template<class _in>
struct rasm_instr0:
    public _in
{
    virtual void to_string(string& str) const override { print(str, get_name(), _t(";")); }
    virtual bool compare(const rasm_instr* that) const override { return strtool::compare(get_name(), that->get_name()) == 0; }
    inline void print(string& str, const gchar* p, const gchar* a) const
    {
        if(p) { str += p; }
        if(a) { str += a; }
    }
    inline void print(string& str, int n, const gchar* a) const
    {
        string s;
        s.format(_t("%08x"), n);
        str += s;
        if(a) { str += a; }
    }
    inline void print(string& str, reg n, const gchar* a) const
    {
        string s;
        s.format(_t("[%s]"), nameofreg(n));
        str += s;
        if(a) { str += a; }
    }
    inline void print(string& str, _bias8_ n, const gchar* a) const
    {
        string s;
        s.format(_t("[s8:%02x]"), (int)n.inner_data());
        str += s;
        if(a) { str += a; }
    }
    inline void print(string& str, _bias16_ n, const gchar* a) const
    {
        string s;
        s.format(_t("[s16:%04x]"), (int)n.inner_data());
        str += s;
        if(a) { str += a; }
    }
    inline void print(string& str, _bias_ n, const gchar* a) const
    {
        string s;
        s.format(_t("[s:%08x]"), (int)n.inner_data());
        str += s;
        if(a) { str += a; }
    }
    inline void print(string& str, _mem_ n, const gchar* a) const
    {
        string s;
        s.format(_t("[m:%08x]"), (int)n.inner_data());
        str += s;
        if(a) { str += a; }
    }
    inline void print(string& str, _reg_ n, const gchar* a) const
    {
        string s;
        s.format(_t("[%s]"), nameofreg(n.inner_data()));
        str += s;
        if(a) { str += a; }
    }
};

template<class c1>
struct rasm_instr1:
    public rasm_instr0<rasm_instr>
{
    typedef c1  type1;
    typedef rasm_instr1<c1> myref;

    c1              _data1;

    rasm_instr1(c1 a) { _data1 = a; }
    virtual void to_string(string& str) const override
    {
        print(str, get_name(), _t("  "));
        print(str, _data1, _t(";"));
    }
    virtual bool compare(const rasm_instr* that) const override
    {
        if(strtool::compare(get_name(), that->get_name()) != 0)
            return false;
        const myref* ptr = static_cast<const myref*>(that);
        return _data1 == ptr->_data1;
    }
};

struct __gs_novtable rasm_link abstract:
    public rasm_instr
{
    virtual bool is_linker() const = 0;
    virtual const string& get_label() const = 0;
};

struct rasm_linker;

struct rasm_label:
    public rasm_instr0<rasm_link>
{
    typedef rasm_linker linker;
    typedef gs::vector<linker*> linker_list;

    string          _label;
    linker_list     _linkers;

public:
    rasm_label(lab l) { _label.assign(l); }
    virtual const gchar* get_name() const override { return _t("label"); }
    virtual bool finalize(vessel& vsl) override { return true; }
    virtual void to_string(string& str) const override;
    virtual bool compare(const rasm_instr* that) const override;
    virtual bool is_linker() const override { return false; }
    virtual const string& get_label() const override { return _label; }

public:
    int get_linker_ctr() const { return (int)_linkers.size(); }
    rasm_linker* get_linker(int i) { return _linkers.at(i); }
    void add_linker(linker* p) { _linkers.push_back(p); }
};

struct rasm_linker:
    public rasm_instr0<rasm_link>
{
    typedef lab type1;
    typedef rasm_linker linker;
    typedef gs::vector<linker*> linker_list;

    string          _label;
    int             _linkeeid;
    linker_list     _affects;

public:
    rasm_linker(lab l): _linkeeid(0) { _label.assign(l); }
    virtual void to_string(string& str) const override;
    virtual bool compare(const rasm_instr* that) const override;
    virtual bool is_linker() const override { return true; }
    virtual const string& get_label() const override { return _label; }

public:
    void set_linkee_index(int i) { _linkeeid = i; }
    int get_linkee_index() const { return _linkeeid; }
    void add_affects(linker* p) { _affects.push_back(p); }
};

template<class c1, class c2>
struct rasm_instr2:
    public rasm_instr1<c1>
{
    typedef c2  type2;
    typedef rasm_instr2<c1, c2> myref;

    c2              _data2;

    rasm_instr2(c1 a, c2 b): rasm_instr1(a) { _data2 = b; }
    virtual void to_string(string& str) const override
    {
        print(str, get_name(), _t("  "));
        print(str, _data1, _t(","));
        print(str, _data2, _t(";"));
    }
    virtual bool compare(const rasm_instr* that) const override
    {
        if(strtool::compare(get_name(), that->get_name()) != 0)
            return false;
        const myref* ptr = static_cast<const myref*>(that);
        return _data1 == ptr->_data1 && _data2 == ptr->_data2;
    }
};

template<class c1, class c2, class c3>
struct rasm_instr3:
    public rasm_instr2<c1, c2>
{
    typedef c3  type3;
    typedef rasm_instr3<c1, c2, c3> myref;

    c3              _data3;

    rasm_instr3(c1 a, c2 b, c3 c): rasm_instr2(a, b) { _data3 = c; }
    virtual void to_string(string& str) const override
    {
        print(str, get_name(), _t("  "));
        print(str, _data1, _t(","));
        print(str, _data2, _t(","));
        print(str, _data3, _t(";"));
    }
    virtual bool compare(const rasm_instr* that) const override
    {
        if(strtool::compare(get_name(), that->get_name()) != 0)
            return false;
        const myref* ptr = static_cast<const myref*>(that);
        return _data1 == ptr->_data1 && _data2 == ptr->_data2 && _data3 == ptr->_data3;
    }
};

#define rasm_reginstr(in) \
    static rasm_instr* make_##in(); \
    void in() { _list.push_back(make_##in()); }

#define rasm_reginstr_i_begin(in, c1) \
    static rasm_instr* make_##in(c1); \
    void in(c1 a) {

#define rasm_reginstr_i_end(in) \
    _list.push_back(make_##in(a)); }

#define rasm_reginstr_i(in, c1) \
    rasm_reginstr_i_begin(in, c1) \
    rasm_reginstr_i_end(in)

#define rasm_reginstr_ii_begin(in, c1, c2) \
    static rasm_instr* make_##in(c1, c2); \
    void in(c1 a, c2 b) {

#define rasm_reginstr_ii_end(in) \
    _list.push_back(make_##in(a, b)); }

#define rasm_reginstr_ii(in, c1, c2) \
    rasm_reginstr_ii_begin(in, c1, c2) \
    rasm_reginstr_ii_end(in)

#define rasm_reginstr_iii_begin(in, c1, c2, c3) \
    static rasm_instr* make_##in(c1, c2, c3); \
    void in(c1 a, c2 b, c3 c) {

#define rasm_reginstr_iii_end(in) \
    _list.push_back(make_##in(a, b, c)); }

#define rasm_reginstr_iii(in, c1, c2, c3) \
    rasm_reginstr_iii_begin(in, c1, c2, c3) \
    rasm_reginstr_iii_end(in)

#define rasm_reginstr_lab(in, c1) \
    static rasm_instr* make_##in(c1); \
    void in(c1 a) \
    { \
        rasm_label* ptr = static_cast<rasm_label*>(make_##in(a)); \
        assert(ptr && "alloc failed."); \
        _list.push_back(ptr); \
        verify(_label.insert(std::make_pair(ptr->get_label(), get_last_instruction_index())).second && "insert failed."); \
    }

#define rasm_reginstr_lnk(in, c1) \
    static rasm_instr* make_##in(c1); \
    void in(c1 a) \
    { \
        rasm_linker* ptr = static_cast<rasm_linker*>(make_##in(a)); \
        assert(ptr && "alloc failed."); \
        _list.push_back(ptr); \
        _link.push_back(ptr); \
    }

/*
 * Only the 8 - bit offsets need this optimization, but DONOT apply such kind of optimizations to immediate operands,
 * that might cause problems like "mov eax, 1" become "mov ah, 1", which was incorrect
 */
#define rasm_optimize_b8(calling, bparam) \
    if((int)bparam.inner_data() > (int)-0x80 && (int)bparam.inner_data() < (int)0x80) \
        return calling;

#define rasm_optimize_b16(calling, bparam) \
    if((int)bparam.inner_data() > (int)-0x8000 && (int)bparam.inner_data() < (int)0x8000) \
        return calling;

/*
 * If a instruction was operated on the eax register, it could be optimized like:
 */
#define rasm_optimize_acc(calling, rparam) \
    if(rparam == _eax) \
        return calling;

/*
 * Optimize for inc or dec
 */
#define rasm_optimize_i1(calling, iparam) \
    if(iparam == 1) \
        return calling;

/*
 * Optimize for immediate operants less than 1 byte. Instructions like loop need this optimization
 */
#define rasm_optimize_i8(calling, iparam) \
    if((int)iparam > (int)-0x80 && (int)iparam < (int)0x80) \
        return calling;

class rasm
{
public:
    typedef gs::vector<rasm_instr*> instr_list;
    typedef gs::vector<rasm_link*> link_list;
    typedef gs::unordered_map<gs::string, int> label_map;       /* label name => index in instruction list */

public:
    rasm() {}
    ~rasm();
    void print(string& str) const;
    int finalize(vessel& vsl);
    int get_last_instruction_index() const;
    //link get_label(lab l) { return _label.find(l); }
    //link no_label() { return _label.end(); }
    //const_link get_label(lab l) const { return _label.find(l); }
    //const_link no_label() const { return _label.end(); }
    //iterator get_current() { return _list.end(); }
    //const_iterator get_current() const { return _list.end(); }
    //iterator get_head() { return _list.begin(); }
    //const_iterator get_head() const { return _list.begin(); }
    //void attach(const_iterator pos, rasm& other);
    //iterator tick() { return _list.size() ? -- _list.end() : _list.end(); }
    //void tick(iterator& i) { i = tick(); }

protected:
    instr_list      _list;
    label_map       _label;
    link_list       _link;

public:
    rasm_reginstr_ii(mov, reg, int8);
    rasm_reginstr_ii(mov, reg, int);
    rasm_reginstr_ii(mov, reg, reg);
    rasm_reginstr_ii(mov, reg, _mem_);
    rasm_reginstr_ii(mov, reg, _bias8_);
    rasm_reginstr_ii_begin(mov, reg, _bias_)
        rasm_optimize_b8(mov(a, (_bias8_)b.inner_data()), b)
        rasm_reginstr_ii_end(mov);
    rasm_reginstr_ii(mov, _mem_, reg);
    rasm_reginstr_ii(mov, _bias8_, reg);
    rasm_reginstr_ii_begin(mov, _bias_, reg)
        rasm_optimize_b8(mov((_bias8_)a.inner_data(), b), a)
        rasm_reginstr_ii_end(mov);
    rasm_reginstr_ii(mov, _mem_, int);
    rasm_reginstr_ii(mov, _bias8_, int);
    rasm_reginstr_ii_begin(mov, _bias_, int)
        rasm_optimize_b8(mov((_bias8_)a.inner_data(), b), a)
        rasm_reginstr_ii_end(mov);
    rasm_reginstr_ii(mov, reg, _reg_);
    rasm_reginstr_ii(mov, _reg_, reg);
    rasm_reginstr_ii(lea, reg, _bias8_);
    rasm_reginstr_ii_begin(lea, reg, _bias_)
        rasm_optimize_b8(lea(a, (_bias8_)b.inner_data()), b)
        rasm_reginstr_ii_end(lea);
    rasm_reginstr(nop);
    rasm_reginstr(hlt);
    rasm_reginstr_i(call, int);
    rasm_reginstr_i(call, reg);
    rasm_reginstr_i(call, _bias_);
    rasm_reginstr_lnk(call, lab);
    rasm_reginstr(ret);
    rasm_reginstr_i(ret, int);
    rasm_reginstr_i(push, reg);
    rasm_reginstr_i(push, int8);
    rasm_reginstr_i(push, int);
    rasm_reginstr_i(push, _mem_);
    rasm_reginstr_i(push, _bias8_);
    rasm_reginstr_i_begin(push, _bias_)
        rasm_optimize_b8(push((_bias8_)a.inner_data()), a)
        rasm_reginstr_i_end(push);
    rasm_reginstr_i(push, _reg_);
    rasm_reginstr_i(pushw, _bias8_);
    rasm_reginstr_i_begin(pushw, _bias_)
        rasm_optimize_b8(pushw((_bias8_)a.inner_data()), a)
        rasm_reginstr_i_end(pushw);
    rasm_reginstr(pushad);
    rasm_reginstr_i(pop, reg);
    rasm_reginstr_i(pop, _mem_);
    rasm_reginstr_i(pop, _bias8_);
    rasm_reginstr_i_begin(pop, _bias_)
        rasm_optimize_b8(pop((_bias8_)a.inner_data()), a)
        rasm_reginstr_i_end(pop);
    rasm_reginstr(popad);
    rasm_reginstr_lab(label, lab);
    rasm_reginstr_i(jmp, _mem_);
    rasm_reginstr_i(jmp, _bias8_);
    rasm_reginstr_i_begin(jmp, _bias_)
        rasm_optimize_b8(jmp((_bias8_)a.inner_data()), a)
        rasm_reginstr_i_end(jmp);
    rasm_reginstr_i(jmp, reg);
    rasm_reginstr_i(jmps, _bias8_);         /* jmp bias from ss */
    rasm_reginstr_i_begin(jmps, _bias_)
        rasm_optimize_b8(jmps((_bias8_)a.inner_data()), a)
        rasm_reginstr_i_end(jmps);
    //rasm_reginstr_i(jmp, lab);
    rasm_reginstr_lnk(jmp, lab);
    rasm_reginstr_i(jes, _bias8_);
    rasm_reginstr_i_begin(jes, _bias_)
        rasm_optimize_b8(jes((_bias8_)a.inner_data()), a)
        rasm_reginstr_i_end(jes);
    rasm_reginstr_i(jnes, _bias8_);
    rasm_reginstr_i_begin(jnes, _bias_)
        rasm_optimize_b8(jnes((_bias8_)a.inner_data()), a)
        rasm_reginstr_i_end(jnes);
    rasm_reginstr_i(loop, int8);
    rasm_reginstr_i_begin(loop, int)        /* I'm not sure about this, can the loop instruction operate on a short? */
        rasm_optimize_i8(loop((int8)a), (a - 2))
        rasm_reginstr_i_end(loop);
    rasm_reginstr_lnk(loop, lab);
    rasm_reginstr(cld);
    rasm_reginstr(std);
    rasm_reginstr(movsb);
    rasm_reginstr(movsw);
    rasm_reginstr(movs);
    rasm_reginstr(rep);
    rasm_reginstr_ii(cmp, _bias_, int);
    rasm_reginstr_ii(cmp, reg, int8);
    rasm_reginstr_ii(cmp, reg, int);
    rasm_reginstr_ii(cmp, reg, _bias8_);
    rasm_reginstr_ii_begin(cmp, reg, _bias_)
        rasm_optimize_b8(cmp(a, (_bias8_)b.inner_data()), b)
        rasm_reginstr_ii_end(cmp);
    rasm_reginstr_ii(cmp, _reg_, int8);
    rasm_reginstr_ii(cmp, _reg_, int);
    rasm_reginstr_ii(cmp, reg, _reg_);
    rasm_reginstr_i(inc, reg);
    rasm_reginstr_i(inc, _mem_);
    rasm_reginstr_i(inc, _bias8_);
    rasm_reginstr_i_begin(inc, _bias_)
        rasm_optimize_b8(inc((_bias8_)a.inner_data()), a)
        rasm_reginstr_i_end(inc);
    rasm_reginstr_ii(add, reg, reg);
    rasm_reginstr_ii(add, reg, int8);
    rasm_reginstr_ii(add, reg, int);
    rasm_reginstr_i(add, int);
    rasm_reginstr_ii(add, reg, _mem_);
    rasm_reginstr_ii(add, reg, _bias8_);
    rasm_reginstr_ii(add, reg, _bias_);
    rasm_reginstr_ii(add, _mem_, int8);
    rasm_reginstr_ii(add, _mem_, int);
    rasm_reginstr_ii(add, _bias8_, int8);
    rasm_reginstr_ii(add, _bias_, int8);
    rasm_reginstr_ii(add, _bias8_, int);
    rasm_reginstr_ii(add, _bias_, int);
    rasm_reginstr_ii(add, _mem_, reg);
    rasm_reginstr_ii(add, _bias8_, reg);
    rasm_reginstr_ii(add, _bias_, reg);
    rasm_reginstr_ii(add, reg, _reg_);
    rasm_reginstr_ii(add, _reg_, reg);
    rasm_reginstr_i(dec, reg);
    rasm_reginstr_i(dec, _mem_);
    rasm_reginstr_i(dec, _bias_);
    rasm_reginstr_ii(sub, reg, reg);
    rasm_reginstr_ii(sub, reg, int);
    rasm_reginstr_i(sub, int);
    rasm_reginstr_ii(sub, reg, _mem_);
    rasm_reginstr_ii(sub, reg, _bias_);
    rasm_reginstr_ii(sub, _mem_, int);
    rasm_reginstr_ii(sub, _bias_, int);
    rasm_reginstr_ii(sub, _mem_, reg);
    rasm_reginstr_ii(sub, _bias_, reg);
    rasm_reginstr_ii(sub, reg, _reg_);
    rasm_reginstr_ii(sub, _reg_, reg);
    rasm_reginstr_iii(imul, reg, reg, int);
    rasm_reginstr_ii(imul, reg, int);
    rasm_reginstr_i(imul, reg);
    rasm_reginstr_ii(imul, reg, _mem_);
    rasm_reginstr_ii(imul, reg, _bias_);
    rasm_reginstr_iii(imul, reg, _mem_, int);
    rasm_reginstr_iii(imul, reg, _bias_, int);
    rasm_reginstr_ii(imul, reg, reg);
    rasm_reginstr_i(imul, _mem_);
    rasm_reginstr_i(imul, _bias_);
    rasm_reginstr_ii(imul, reg, _reg_);
    rasm_reginstr_i(idiv, reg);
    rasm_reginstr_i(idiv, _mem_);
    rasm_reginstr_i(idiv, _bias_);
    rasm_reginstr_i(idiv, _reg_);
    rasm_reginstr(cdq);
    rasm_reginstr_i(neg, reg);
    rasm_reginstr_i(not, reg);
    rasm_reginstr_ii(xor, reg, reg);
    rasm_reginstr_ii(xor, reg, int);
    rasm_reginstr_ii(xor, reg, _bias_);
    rasm_reginstr_ii(xor, reg, _reg_);
    rasm_reginstr_i(shl, reg);
    rasm_reginstr_ii(shl, reg, int);
    rasm_reginstr_i(shr, reg);
    rasm_reginstr_ii(shr, reg, int);
    //rasm_reginstr
    //setcc
    rasm_reginstr_ii(and, reg, int);
    rasm_reginstr_ii(and, reg, _bias_);
    rasm_reginstr_ii(and, reg, _reg_);
    rasm_reginstr_ii(or , reg, int);
    rasm_reginstr_ii(or , reg, _bias_);
    rasm_reginstr_ii(or , reg, _reg_);
};

#endif

__rathen_end__

#endif

ext: h

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
/*
 * Copyright 2016 Maikel Nadolski
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0

 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#ifndef HMM_ALGORITHM_FORWARD_H_
#define HMM_ALGORITHM_FORWARD_H_

#include <utility>
#include <boost/iterator/iterator_facade.hpp>
#include <gsl_assert.h>

#include "maikel/hmm/hidden_markov_model.h"

namespace maikel { namespace hmm {

  template <class InputIter, class T>
    class forward_range_fn {
      public:
        using model       = hidden_markov_model<T>;
        using row_vector  = typename model::row_vector;
        using matrix      = typename model::matrix;
        using symbol_type = typename std::iterator_traits<InputIter>::value_type;

        forward_range_fn() = delete;

        forward_range_fn(InputIter seq_it, InputIter seq_end, model const& hmm)
        :  hmm_{&hmm}, seq_it_{seq_it}, seq_end_{seq_end},
          alpha_{0,row_vector(hmm.states())}, prev_alpha_(hmm.states())
        {
          if (seq_it != seq_end)
            initial_coefficients(*seq_it);
        }

        class iterator
        : public boost::iterator_facade<
              iterator, std::pair<T, row_vector>, std::input_iterator_tag,
              std::pair<T, row_vector> const&
         > {
          public:
            iterator() = default;
          private:
            friend class forward_range_fn;
            friend class boost::iterator_core_access;

            iterator(forward_range_fn& parent)
            : parent_{parent ? &parent : nullptr} {}

            forward_range_fn* parent_ = nullptr;

            std::pair<T, row_vector> const& dereference() const
            {
              Expects(parent_ && *parent_);
              return parent_->alpha_;
            }

            void increment()
            {
              Expects(parent_ && *parent_);
              if (!parent_->next())
                parent_ = nullptr;
            }

            bool equal(iterator other) const noexcept
            {
              return parent_ == other.parent_;
            }
        };

        operator bool() const noexcept
        {
          return seq_it_ != seq_end_;
        }

        iterator begin() noexcept
        {
          return {*this};
        }

        iterator end() noexcept
        {
          return {};
        }

      private:
        model const* hmm_; // not owning
        InputIter seq_it_, seq_end_;
        std::pair<T, row_vector> alpha_;
        row_vector prev_alpha_;

        void initial_coefficients(symbol_type s)
        {
          using size_type = typename hidden_markov_model<T>::size_type;
          matrix const& B  = hmm_->symbol_probabilities();
          row_vector const& pi = hmm_->initial_distribution();

          // check pre conditions
          size_type states = pi.size();
          size_type ob = gsl::narrow<size_type>(s);
          Expects(B.rows() == states);
          Expects(0 <= ob && ob < B.cols());
          Expects(alpha_.second.size() == states);

          // initial formula
          T& scaling = alpha_.first;
          row_vector& alpha = alpha_.second;
          scaling = 0.0;
          for (size_type i = 0; i < states; ++i) {
            alpha(i) = pi(i)*B(i,ob);
            scaling += alpha(i);
          }
          scaling = scaling ? 1/scaling : 0;
          alpha *= scaling;

          // check post conditions
          Ensures((!scaling && almost_equal<T>(alpha.sum(), 0.0)) ||
                  ( scaling && almost_equal<T>(alpha.sum(), 1.0))    );
        }

        void recursion_advance(symbol_type s)
        {
          using size_type = typename hidden_markov_model<T>::size_type;
          matrix const& A = hmm_->transition_matrix();
          matrix const& B = hmm_->symbol_probabilities();
          prev_alpha_.swap(alpha_.second);

          // check pre conditions
          size_type states = A.rows();
          Expects(A.cols() == states);
          Expects(B.rows() == states);
          Expects(prev_alpha_.size() == states);
          size_type ob = gsl::narrow<size_type>(s);
          Expects(0 <= ob && ob < B.cols());

          // recursion formula
          T& scaling = alpha_.first;
          row_vector& alpha = alpha_.second;
          scaling = 0.0;
          for (size_type j = 0; j < states; ++j) {
            alpha(j) = 0.0;
            for (size_type i = 0; i < states; ++i)
              alpha(j) += prev_alpha_(i)*A(i,j);
            alpha(j) *= B(j, ob);
            scaling += alpha(j);
          }
          scaling = scaling ? 1/scaling : 0;
          alpha *= scaling;

          // post conditions
          Ensures((!scaling && almost_equal<T>(alpha.sum(), 0.0)) ||
                  ( scaling && almost_equal<T>(alpha.sum(), 1.0))    );
        }

        bool next()
        {
          Expects(*this);
          ++seq_it_;
          if (seq_it_ == seq_end_)
            return false;
          recursion_advance(*seq_it_);
          return true;
        }
    };

  template <class InputIter, class T>
    forward_range_fn<InputIter, T>
  forward(InputIter begin, InputIter end, hidden_markov_model<T> const& hmm)
  {
    return {begin, end, hmm};
  }

} // namespace hmm
} // namespace maikel

#endif /* HMM_ALGORITHM_FORWARD_H_ */

max_stars_repo_name: maikel/hidden-markov-model

max_issues_repo_name: maikel/Hidden-Markov-Model

max_forks_repo_name: maikel/Hidden-Markov-Model


text:
#pragma once

#include "common.h"
#include "multibase.h"
#include "uvarint.h"
#include <gsl/gsl>
#include <map>

namespace multiformats {

    enum hash_t {
        dynamic_hash = -1,
        sha1,
        sha2_256
    };

    namespace details {

        struct hashimpl {
            hash_t      key;
            const char* name;
            uint32_t    code;
            int32_t     len;
        };

        constexpr hashimpl _HashTable[] = {
            { dynamic_hash, "dynamic_hash", 0x00,  0 },
            { sha1,         "sha1",         0x11, 20 },
            { sha2_256,     "sha2_256",     0x12, 32 },
        };

        constexpr int find_hashimpl_by_key(hash_t code) {
            for (auto i = 0; i < _countof(_HashTable); i++)
                if (_HashTable[i].key == code) return i;
            return 0;
        }

        constexpr int find_hashimpl_by_code(uint32_t code) {
            for (auto i = 0; i < _countof(_HashTable); i++)
                if (_HashTable[i].code == code) return i;
            return 0;
        }

        template <hash_t _Hash = dynamic_hash, int _Index = find_hashimpl_by_key(_Hash)>
        class hashcode_type
        {
        public:
            static_assert(_Hash != dynamic_hash, "dynamic_hash is not allowed here");
            static_assert(_Index > 0, "The hash_t is not implemented");

            constexpr hashcode_type(hash_t hash) { Expects(hash == _Hash); }

            constexpr hash_t      hash()  const { return _HashTable[_Index].key; }
            constexpr const char* name()  const { return _HashTable[_Index].name; }
            constexpr uint32_t    code()  const { return _HashTable[_Index].code; }
            constexpr int32_t     len()   const { return _HashTable[_Index].len; }
            constexpr int         index() const { return _Index; }
        };

        template <>
        class hashcode_type<dynamic_hash>
        {
        public:
            explicit constexpr hashcode_type(hash_t hash) : _index(find_hashimpl_by_key(hash)) { Expects(_index > 0); }

            const hash_t      hash()  const { return _HashTable[_index].key; }
            const char*       name()  const { return _HashTable[_index].name; }
            const uint32_t    code()  const { return _HashTable[_index].code; }
            const int32_t     len()   const { return _HashTable[_index].len; }
            const int         index() const { return _index; }
        private:
            const int _index;
        };

    }

    // A strongly typed byte buffer that contains an immutable representation of the output digest for a specific hash algorithm
    // known at compile-time.
    // The specialization with dynamic_hash is used when the hash function is not known at compile-time.

    template <hash_t _Hash = dynamic_hash>
    class digest_buffer
    {
    public:
        // Construct empty
        digest_buffer() : _hash(_Hash) {}
        digest_buffer(hash_t hash) : _hash(hash) {}

        // Construct by copying _Right
        //    - from bufferview_t
        digest_buffer(bufferview_t _Right) : _buffer(std::begin(_Right), std::end(_Right)), _hash(_Hash)
        {
            static_assert(_Hash != dynamic_hash, "dynamic_hash is not allowed here");
            Expects(!_Right.empty());
            Expects(_hash.len() < 0 || _Right.size() == _hash.len());
        }
        digest_buffer(hash_t hash, bufferview_t _Right) : _buffer(std::begin(_Right), std::end(_Right)), _hash(hash)
        {
            Expects(!_Right.empty());
            Expects(_hash.len() < 0 || _Right.size() == _hash.len());
        }

        //    - from digest_buffer<>
        digest_buffer(const digest_buffer<_Hash>& _Right) : _buffer(_Right._buffer), _hash(_Right.hash())
        {
            Expects(!_Right.empty());
        }

        template <hash_t _RightHash>
        digest_buffer(const digest_buffer<_RightHash>& _Right) : _buffer(_Right._buffer), _hash(_Right.hash())
        {
            Expects(!_Right.empty());
        }

        // Construct by moving _Right
        digest_buffer(buffer_t&& _Right) : _buffer(std::move(_Right)), _hash(_Hash)
        {
            static_assert(_Hash != dynamic_hash, "dynamic_hash is not allowed here");
        }
        digest_buffer(hash_t hash, buffer_t&& _Right) : _buffer(std::move(_Right)), _hash(hash) {}

        digest_buffer(digest_buffer<_Hash>&& _Right) : _buffer(std::move(_Right._buffer)), _hash(_Right.hash()) {}
        template <hash_t _RightHash>
        digest_buffer(digest_buffer<_RightHash>&& _Right) : _string(std::move(_Right._buffer)), _hash(_Right.hash()) {}

        // Assign by copying _Right
        // Note: an digest_buffer<> is immutable once initialized
        digest_buffer<_Hash>& operator =(bufferview_t _Right)
        {
            Expects(_buffer.empty() && !_Right.empty());
            _buffer = _Right;
            return *this;
        }

        digest_buffer<_Hash>& operator =(const digest_buffer<_Hash>& _Right)
        {
            Expects(_Right.hash() == hash());
            Expects(_buffer.empty() && !_Right.empty());
            _buffer = _Right._buffer;
            return *this;
        }

        template <hash_t _RightHash>
        digest_buffer<_Hash>& operator =(const digest_buffer<_RightHash>& _Right)
        {
            static_assert(_RightHash == _Hash || _RightHash == dynamic_hash || _Hash == dynamic_hash, "Mismatch between codes.");
            Expects(_Right.hash() == hash());
            Expects(_buffer.empty() && !_Right.empty());
            _buffer = _Right._buffer;
            return *this;
        }

        // Assign by moving _Right
        // Note: an digest_buffer<> is immutable once initialized and cannot be emptied
        digest_buffer<_Hash>& operator =(bufferview_t&& _Right)
        {
            Expects(_buffer.empty() && !_Right.empty());
            _buffer = std::move(_Right);
            return *this;
        }

        digest_buffer<_Hash>& operator =(digest_buffer<_Hash>&& _Right)
        {
            Expects(_Right.hash() == hash());
            Expects(_buffer.empty() && !_Right.empty());
            _buffer = std::move(_Right._buffer);
            return *this;
        }

        template <hash_t _RightHash>
        digest_buffer<_Hash>& operator =(digest_buffer<_RightHash>&& _Right)
        {
            static_assert(_RightHash == _Hash || _RightHash == dynamic_hash || _Hash == dynamic_hash, "Mismatch between codes.");
            Expects(_Right.hash() == hash());
            Expects(_buffer.empty() && !_Right.empty());
            _buffer = std::move(_Right._buffer);
            return *this;
        }

        auto hash()  const { return _hash.hash(); }
        auto data()  const { return _buffer; }
        auto code()  const { return _hash.code(); }

        auto empty() const { return _buffer.empty(); }
        auto size()  const { return _buffer.size(); }
        auto begin() const { return _buffer.begin(); }
        auto end()   const { return _buffer.end(); }

    private:
        const details::hashcode_type<_Hash> _hash;
        buffer_t _buffer;

        template <hash_t _FriendHash> friend class digest_buffer;
    };

    // operator==
    template <hash_t _HashLeft, hash_t _HashRight>
    bool operator==(const digest_buffer<_HashLeft>& _Left, const digest_buffer<_HashRight>& _Right)
    {
        return (_Left.hash() == _Right.hash()) && (_Left.data() == _Right.data());
    }

    // A string that contains an encoded string in base _Base of a digest with algo _Hash
    template <hash_t _Hash = dynamic_hash, base_t _Base = dynamic_base>
    class digest_string : public encoded_string<_Base>
    {
        using string_t = encoded_string<_Base>;

    public:
        // Construct empty
        digest_string() : string_t(), _hash(_Hash) {}
        digest_string(hash_t hash) : string_t(), _hash(hash) {}
        digest_string(base_t base) : string_t(base), _hash(_Hash) {}
        digest_string(hash_t hash, base_t base) : string_t(base), _hash(hash) {}

        // Construct by copying _Right
        //    - from string
        digest_string(stringview_t _Right) : string_t(_Right), _hash(_Hash)
        {
            static_assert(_Hash != dynamic_hash, "dynamic_hash is not allowed here");
            Expects(!_Right.empty());
            Expects(_hash.len() < 0 || decode(this->base(), _Right).size() == _hash.len());
        }
        digest_string(hash_t hash, stringview_t _Right) : string_t(_Right), _hash(hash)
        {
            Expects(!_Right.empty());
            Expects(_hash.len() < 0 || decode(this->base(), _Right).size() == _hash.len());
        }
        digest_string(base_t base, stringview_t _Right) : string_t(base, _Right), _hash(_Hash)
        {
            static_assert(_Hash != dynamic_hash, "dynamic_hash is not allowed here");
            Expects(!_Right.empty());
            Expects(_hash.len() < 0 || decode(this->base(), _Right).size() == _hash.len());
        }
        digest_string(hash_t hash, base_t base, stringview_t _Right) : string_t(base, _Right), _hash(hash)
        {
            Expects(!_Right.empty());
            Expects(_hash.len() < 0 || decode(this->base(), _Right).size() == _hash.len());
        }

        digest_string(const char* _Right) : digest_string(gsl::ensure_z(_Right)) {}
        digest_string(hash_t hash, const char* _Right) : digest_string(hash, gsl::ensure_z(_Right)) {}
        digest_string(base_t base, const char* _Right) : digest_string(base, gsl::ensure_z(_Right)) {}
        digest_string(hash_t hash, base_t base, const char* _Right) : digest_string(hash, base, gsl::ensure_z(_Right)) {}

        //    - from encoded_string<>
        digest_string(const encoded_string<_Base>& _Right) : string_t(_Right), _hash(_Hash)
        {
            static_assert(_Hash != dynamic_hash, "dynamic_hash is not allowed here");
            Expects(!_Right.empty());
            Expects(_hash.len() < 0 || decode(_Right).size() == _hash.len());
        }
        digest_string(hash_t hash, const encoded_string<_Base>& _Right) : string_t(_Right), _hash(hash)
        {
            Expects(!_Right.empty());
            Expects(_hash.len() < 0 || decode(_Right).size() == _hash.len());
        }

        template <base_t _RightBase>
        digest_string(const encoded_string<_RightBase>& _Right) : string_t(_Right), _hash(_Hash)
        {
            static_assert(_RightBase == _Base || _RightBase == dynamic_base || _Base == dynamic_base, "Mismatch between codes.");
            static_assert(_Hash != dynamic_hash, "dynamic_hash is not allowed here");
            Expects(!_Right.empty());
            Expects(_hash.len() < 0 || decode(_Right).size() == _hash.len());
        }
        template <base_t _RightBase>
        digest_string(hash_t hash, const encoded_string<_RightBase>& _Right) : string_t(_Right), _hash(hash)
        {
            static_assert(_RightBase == _Base || _RightBase == dynamic_base || _Base == dynamic_base, "Mismatch between codes.");
            Expects(!_Right.empty());
            Expects(_hash.len() < 0 || decode(_Right).size() == _hash.len());
        }

        //    - from digest_string<>
        digest_string(const digest_string<_Hash, _Base>& _Right) : string_t(_Right), _hash(_Right.hash())
        {
            Expects(!_Right.empty());
        }
        template <hash_t _RightHash>
        digest_string(const digest_string<_RightHash, _Base>& _Right) : string_t(_Right), _hash(_Right.hash())
        {
            static_assert(_RightHash == _Hash || _RightHash == dynamic_hash || _Hash == dynamic_hash, "Mismatch between codes.");
            Expects(!_Right.empty());
        }
        template <base_t _RightBase>
        digest_string(const digest_string<_Hash, _RightBase>& _Right) : string_t(_Right), _hash(_Right.hash())
        {
            static_assert(_RightBase == _Base || _RightBase == dynamic_base || _Base == dynamic_base, "Mismatch between codes.");
            Expects(!_Right.empty());
        }
        template <hash_t _RightHash, base_t _RightBase>
        digest_string(const digest_string<_RightHash, _RightBase>& _Right) : string_t(_Right), _hash(_Right.hash())
        {
            static_assert(_RightHash == _Hash || _RightHash == dynamic_hash || _Hash == dynamic_hash, "Mismatch between codes.");
            static_assert(_RightBase == _Base || _RightBase == dynamic_base || _Base == dynamic_base, "Mismatch between codes.");
            Expects(!_Right.empty());
        }

        // Assign by copying _Right
        // Note: an digest_string<> is immutable once initialized
        digest_string<_Hash, _Base>& operator =(stringview_t _Right)
        {
            Expects(_hash.len() < 0 || decode(this->base(), _Right).size() == _hash.len());
            *((string_t*)this) = _Right;
            return *this;
        }
        digest_string<_Hash, _Base>& operator =(const char* _Right)
        {
            return operator=(gsl::ensure_z(_Right));
        }

        digest_string<_Hash, _Base>& operator =(const encoded_string<_Base>& _Right)
        {
            Expects(_hash.len() < 0 || decode(_Right).size() == _hash.len());
            *((string_t*)this) = _Right;
            return *this;
        }
        template <base_t _RightBase>
        digest_string<_Hash, _Base>& operator =(const encoded_string<_RightBase>& _Right)
        {
            static_assert(_RightBase == _Base || _RightBase == dynamic_base || _Base == dynamic_base, "Mismatch between codes.");
            Expects(_hash.len() < 0 || decode(_Right).size() == _hash.len());
            *((string_t*)this) = _Right;
            return *this;
        }

        digest_string<_Hash, _Base>& operator =(const digest_string<_Hash, _Base>& _Right)
        {
            Expects(_Right.hash() == hash());
            *this = _Right.str();
            return *this;
        }
        template <hash_t _RightHash, base_t _RightBase>
        digest_string<_Hash, _Base>& operator =(const digest_string<_RightHash, _RightBase>& _Right)
        {
            static_assert(_RightHash == _Hash || _RightHash == dynamic_hash || _Hash == dynamic_hash, "Mismatch between codes.");
            static_assert(_RightBase == _Base || _RightBase == dynamic_base || _Base == dynamic_base, "Mismatch between codes.");
            Expects(_Right.hash() == hash());
            *this = _Right.str();
            return *this;
        }

        hash_t hash() const { return _hash.hash(); }

    private:
        const details::hashcode_type<_Hash> _hash;
    };

    //
    // Encode a digest into a _Base encoded_string
    //
    template <hash_t _Hash, base_t _Base>
    digest_string<_Hash, _Base> encode(const digest_buffer<_Hash>& digest)
    {
        Expects(!digest.empty());
        return encode<_Base>(digest);
    }

    template <hash_t _Hash>
    digest_string<_Hash> encode(base_t base, const digest_buffer<_Hash>& digest)
    {
        Expects(!digest.empty());
        return encode(base, digest);
    }

    inline digest_string<> encode(hash_t hash, base_t base, bufferview_t digest)
    {
        Expects(!digest.empty());
        return { hash, encode(base, digest) };
    }

    template <hash_t _Hash, base_t _Base>
    digest_buffer<_Hash> decode(digest_string<_Hash, _Base> digest)
    {
        return decode(static_cast<encoded_string<_Base>>(digest));
    }

    class multihash
    {
    public:
        multihash(bufferview_t mhview)
        {
            auto posIt = mhview.begin();

            // parse and check hash code
            auto code = uint32_t{};
            posIt = uvarint::decode(posIt, mhview.end(), &code);

            auto index = details::find_hashimpl_by_code(code);
            if (index == 0) throw std::invalid_argument("Failed to parse multihash buffer: wrong hash code");
            _hash = details::_HashTable[index].key;

            // parse and check size
            posIt = uvarint::decode(posIt, mhview.end(), &_size);

            if (static_cast<ptrdiff_t>(_size) != mhview.end() - posIt) throw std::invalid_argument("Failed to parse multihash buffer: wrong size");
            if (details::_HashTable[index].len > 0 && _size != details::_HashTable[index].len) throw std::invalid_argument("Failed to parse multihash buffer: wrong size");

            // Seems legit
            _data.assign(mhview.begin(), mhview.end());
        }

        hash_t       hash()   const { return _hash; }
        size_t       size()   const { return _size; }
        bufferview_t digest() const { return bufferview_t{ _data }.last(_size); }

        bufferview_t data()   const { return _data; }

    private:
        hash_t _hash;
        size_t _size;
        buffer_t _data;
    };

    // Create a multihash from a digest_buffer
    template <hash_t _Hash>
    multihash to_multihash(digest_buffer<_Hash> digest) {
        auto mh = buffer_t{};

        uvarint::encode(digest.code(), std::back_inserter(mh));
        uvarint::encode(digest.size(), std::back_inserter(mh));
        mh.insert(mh.end(), digest.begin(), digest.end());

        return { mh };
    }

    // Create a multihash from a digest_string
    template <hash_t _Hash, base_t _Base>
    multihash to_multihash(digest_string<_Hash, _Base> digest) { return to_multihash(decode(digest)); }

    // Create a multihash from a raw digest_string
    template <hash_t _Hash>
    multihash to_multihash(bufferview_t digest) { return to_multihash(digest_buffer<_Hash>{ digest }); }
    inline multihash to_multihash(hash_t hash, bufferview_t digest) { return to_multihash(digest_buffer<>{ hash, digest }); }
}

max_stars_repo_name: cedrou/cpp-multiformats

max_issues_repo_name: cedrou/multiformats

max_forks_repo_name: cedrou/multiformats


text:
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>
#include <cairo.h>
#include <gtk/gtk.h>

double globx;

double analytic(double w0, double t);

double spring_euler(double w0, double x);
void test_euler(GtkWidget *darea);

int func(double t, const double x[], double f[], void *params);
int jac(double t, const double x[], double *dfdy, double dfdt[], void *params);
int test_gsl(GtkWidget *darea);

void do_drawing(cairo_t *cr);
gboolean on_draw_event(GtkWidget *widget, cairo_t *cr, gpointer user_data);

gboolean terminate(GThread *thread);

void *thread_func(GtkWidget *entry) {
  //test_euler(entry);
  test_gsl(entry);

  gdk_threads_add_idle((GSourceFunc)terminate, g_thread_self());
  return NULL;
}

int main(int argc, char *argv[]) {
	GtkWidget *window;
	GtkWidget *darea;

    gtk_init(&argc, &argv);

	window = gtk_window_new(GTK_WINDOW_TOPLEVEL);

	darea = gtk_drawing_area_new();
	gtk_container_add(GTK_CONTAINER(window), darea);

	g_signal_connect(G_OBJECT(darea), "draw", G_CALLBACK(on_draw_event), NULL);
	g_signal_connect(window, "destroy",	G_CALLBACK(gtk_main_quit), NULL);

	gtk_window_set_position(GTK_WINDOW(window), GTK_WIN_POS_CENTER);
	gtk_window_set_default_size(GTK_WINDOW(window), 400, 300);
	gtk_window_set_title(GTK_WINDOW(window), "Spring simulation");

	gtk_widget_show_all(window);

  	g_thread_new("eval", (GThreadFunc)thread_func, darea);

	gtk_main();

	return 0;
}

double analytic(double w0, double t) {
	return cos(w0 * t);
}

// euler

double spring_euler(double w0, double x) {
	return -w0 * w0 * x;
}

void test_euler(GtkWidget *darea) {
	const int n = 1000;

	// x - position
	// v - speed
	// t - time
	// dt - step
	// w0 = sqrt(k/m)
	double x[n], v[n], t[n], dt;

	const double w0 = 1.0;
	x[0] = 1.0; v[0] = 0.0; t[0] = 0.0;
	dt = 0.05;

	int i;
	for (i = 0; i < n - 1 ; i++) {
		x[i+1] = x[i] + v[i] * dt;
		v[i+1] = v[i] + spring_euler(w0, x[i]) * dt;
		t[i+1] = t[i] + dt;
	}

	for (i = 0; i < n; i++) {
		// draw
		globx = x[i];
		gtk_widget_queue_draw(darea);
		printf("t: %f\tx: %f\tanalytic: \t%f\n", t[i], x[i], analytic(w0, t[i]));
		g_usleep(10000);
	}
}

// gsl

int func(double t, const double x[], double f[], void *params) {
	// x' = v
	// v' = -w0^2 * x
	double w0 = *(double *) params;
	f[0] = x[1];
	f[1] = -w0 * w0 * x[0];
	return GSL_SUCCESS;
}

int jac(double t, const double x[], double *dfdy, double dfdt[], void *params) {
	double w0 = *(double *)params;
	gsl_matrix_view dfdy_mat = gsl_matrix_view_array(dfdy, 2, 2);
	gsl_matrix * m = &dfdy_mat.matrix;
	gsl_matrix_set(m, 0, 0, 0.0);
	gsl_matrix_set(m, 0, 1, 1.0);
	gsl_matrix_set(m, 1, 0, -w0 * w0);
	gsl_matrix_set(m, 1, 1, 0.0);
	dfdt[0] = 0.0;
	dfdt[1] = 0.0;
	return GSL_SUCCESS;
}

int test_gsl(GtkWidget *darea) {
	double w0 = 1.0;
	gsl_odeiv2_system sys = { func, jac, 2, &w0 };

	gsl_odeiv2_driver *d = gsl_odeiv2_driver_alloc_y_new(&sys, gsl_odeiv2_step_rk4, 0.05, 1e-6, 1e-6);

	double t = 0.0;
	double x[2] = { 1.0, 0.0 };
	int i, s;

	for (i = 0; i < 1000; i++) {
		s = gsl_odeiv2_driver_apply_fixed_step(d, &t, 0.05, 1, x);

		if (s != GSL_SUCCESS) {
			printf ("error: driver returned %d\n", s);
			break;
		}

		// draw
		globx = x[0];
		gtk_widget_queue_draw(darea);
		printf("t: %.4f\tanalytic: %.8f\tx:%.8f\n", t, analytic(w0, t), x[0]);
		g_usleep(10000);
	}

	gsl_odeiv2_driver_free(d);
	return s;
}

// drawing

gboolean on_draw_event(GtkWidget *widget, cairo_t *cr, gpointer user_data) {
	do_drawing(cr);

	return FALSE;
}

void do_drawing(cairo_t *cr) {
	cairo_set_source_rgb(cr, 0, 0, 0);
	cairo_set_line_width(cr, 1);

	int k;
	int prevy = 0;
	cairo_move_to(cr, 200, 0);
	for (k = 0; k < 20; k++) {
		int j = k % 2 == 0 ? -1 : 1;
		cairo_line_to(cr, 200 + j * 20, k*5*(globx + 1) + prevy);
		prevy = k*2;
	}
	cairo_line_to(cr, 200, 100*(globx + 1) + prevy);
	//cairo_stroke(cr);
	cairo_rectangle(cr, 180, 100*(globx + 1) + prevy, 40, 40);
	//cairo_stroke_preserve(cr);
	//cairo_fill(cr);
	cairo_stroke(cr);
}

// threads

gboolean terminate(GThread *thread) {
  g_thread_join(thread);

  //g_thread_unref(thread);

  return FALSE;
}

ext: c

max_stars_repo_name: Hoobie/spring-simulation

max_issues_repo_name: Hoobie/spring-simulation

max_forks_repo_name: Hoobie/spring-simulation


text:
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_cblas.h>

#include "tests.h"

void
test_spmv (void) {
const double flteps = 1e-4, dbleps = 1e-6;
  {
   int order = 101;
   int uplo = 121;
   float alpha = 0.1f;
   float beta = -0.3f;
   int N = 2;
   float A[] = { -0.174f, 0.878f, 0.478f };
   float X[] = { 0.503f, 0.313f };
   int incX = -1;
   float Y[] = { -0.565f, -0.109f };
   int incY = -1;
   float y_expected[] = { 0.221025f, 0.0714172f };
   cblas_sspmv(order, uplo, N, alpha, A, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(Y[i], y_expected[i], flteps, "sspmv(case 1134)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 121;
   float alpha = 0.1f;
   float beta = -0.3f;
   int N = 2;
   float A[] = { -0.174f, 0.878f, 0.478f };
   float X[] = { 0.503f, 0.313f };
   int incX = -1;
   float Y[] = { -0.565f, -0.109f };
   int incY = -1;
   float y_expected[] = { 0.221025f, 0.0714172f };
   cblas_sspmv(order, uplo, N, alpha, A, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(Y[i], y_expected[i], flteps, "sspmv(case 1135)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 122;
   float alpha = 0.1f;
   float beta = -0.3f;
   int N = 2;
   float A[] = { -0.174f, 0.878f, 0.478f };
   float X[] = { 0.503f, 0.313f };
   int incX = -1;
   float Y[] = { -0.565f, -0.109f };
   int incY = -1;
   float y_expected[] = { 0.221025f, 0.0714172f };
   cblas_sspmv(order, uplo, N, alpha, A, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(Y[i], y_expected[i], flteps, "sspmv(case 1136)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 122;
   float alpha = 0.1f;
   float beta = -0.3f;
   int N = 2;
   float A[] = { -0.174f, 0.878f, 0.478f };
   float X[] = { 0.503f, 0.313f };
   int incX = -1;
   float Y[] = { -0.565f, -0.109f };
   int incY = -1;
   float y_expected[] = { 0.221025f, 0.0714172f };
   cblas_sspmv(order, uplo, N, alpha, A, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(Y[i], y_expected[i], flteps, "sspmv(case 1137)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 121;
   float alpha = 0.1f;
   float beta = -0.3f;
   int N = 2;
   float A[] = { -0.174f, 0.878f, 0.478f };
   float X[] = { 0.503f, 0.313f };
   int incX = -1;
   float Y[] = { -0.565f, -0.109f };
   int incY = -1;
   float y_expected[] = { 0.221025f, 0.0714172f };
   cblas_sspmv(order, uplo, N, alpha, A, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(Y[i], y_expected[i], flteps, "sspmv(case 1138)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 121;
   float alpha = 0.1f;
   float beta = -0.3f;
   int N = 2;
   float A[] = { -0.174f, 0.878f, 0.478f };
   float X[] = { 0.503f, 0.313f };
   int incX = -1;
   float Y[] = { -0.565f, -0.109f };
   int incY = -1;
   float y_expected[] = { 0.221025f, 0.0714172f };
   cblas_sspmv(order, uplo, N, alpha, A, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(Y[i], y_expected[i], flteps, "sspmv(case 1139)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 122;
   float alpha = 0.1f;
   float beta = -0.3f;
   int N = 2;
   float A[] = { -0.174f, 0.878f, 0.478f };
   float X[] = { 0.503f, 0.313f };
   int incX = -1;
   float Y[] = { -0.565f, -0.109f };
   int incY = -1;
   float y_expected[] = { 0.221025f, 0.0714172f };
   cblas_sspmv(order, uplo, N, alpha, A, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(Y[i], y_expected[i], flteps, "sspmv(case 1140)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 122;
   float alpha = 0.1f;
   float beta = -0.3f;
   int N = 2;
   float A[] = { -0.174f, 0.878f, 0.478f };
   float X[] = { 0.503f, 0.313f };
   int incX = -1;
   float Y[] = { -0.565f, -0.109f };
   int incY = -1;
   float y_expected[] = { 0.221025f, 0.0714172f };
   cblas_sspmv(order, uplo, N, alpha, A, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(Y[i], y_expected[i], flteps, "sspmv(case 1141)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 121;
   double alpha = -1;
   double beta = 0.1;
   int N = 2;
   double A[] = { -0.181, -0.071, -0.038 };
   double X[] = { -0.015, 0.132 };
   int incX = -1;
   double Y[] = { -0.449, -0.219 };
   int incY = -1;
   double y_expected[] = { -0.036098, 9.27e-04 };
   cblas_dspmv(order, uplo, N, alpha, A, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(Y[i], y_expected[i], dbleps, "dspmv(case 1142)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 121;
   double alpha = -1;
   double beta = 0.1;
   int N = 2;
   double A[] = { -0.181, -0.071, -0.038 };
   double X[] = { -0.015, 0.132 };
   int incX = -1;
   double Y[] = { -0.449, -0.219 };
   int incY = -1;
   double y_expected[] = { -0.036098, 9.27e-04 };
   cblas_dspmv(order, uplo, N, alpha, A, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(Y[i], y_expected[i], dbleps, "dspmv(case 1143)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 122;
   double alpha = -1;
   double beta = 0.1;
   int N = 2;
   double A[] = { -0.181, -0.071, -0.038 };
   double X[] = { -0.015, 0.132 };
   int incX = -1;
   double Y[] = { -0.449, -0.219 };
   int incY = -1;
   double y_expected[] = { -0.036098, 9.27e-04 };
   cblas_dspmv(order, uplo, N, alpha, A, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(Y[i], y_expected[i], dbleps, "dspmv(case 1144)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 122;
   double alpha = -1;
   double beta = 0.1;
   int N = 2;
   double A[] = { -0.181, -0.071, -0.038 };
   double X[] = { -0.015, 0.132 };
   int incX = -1;
   double Y[] = { -0.449, -0.219 };
   int incY = -1;
   double y_expected[] = { -0.036098, 9.27e-04 };
   cblas_dspmv(order, uplo, N, alpha, A, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(Y[i], y_expected[i], dbleps, "dspmv(case 1145)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 121;
   double alpha = -1;
   double beta = 0.1;
   int N = 2;
   double A[] = { -0.181, -0.071, -0.038 };
   double X[] = { -0.015, 0.132 };
   int incX = -1;
   double Y[] = { -0.449, -0.219 };
   int incY = -1;
   double y_expected[] = { -0.036098, 9.27e-04 };
   cblas_dspmv(order, uplo, N, alpha, A, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(Y[i], y_expected[i], dbleps, "dspmv(case 1146)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 121;
   double alpha = -1;
   double beta = 0.1;
   int N = 2;
   double A[] = { -0.181, -0.071, -0.038 };
   double X[] = { -0.015, 0.132 };
   int incX = -1;
   double Y[] = { -0.449, -0.219 };
   int incY = -1;
   double y_expected[] = { -0.036098, 9.27e-04 };
   cblas_dspmv(order, uplo, N, alpha, A, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(Y[i], y_expected[i], dbleps, "dspmv(case 1147)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 122;
   double alpha = -1;
   double beta = 0.1;
   int N = 2;
   double A[] = { -0.181, -0.071, -0.038 };
   double X[] = { -0.015, 0.132 };
   int incX = -1;
   double Y[] = { -0.449, -0.219 };
   int incY = -1;
   double y_expected[] = { -0.036098, 9.27e-04 };
   cblas_dspmv(order, uplo, N, alpha, A, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(Y[i], y_expected[i], dbleps, "dspmv(case 1148)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 122;
   double alpha = -1;
   double beta = 0.1;
   int N = 2;
   double A[] = { -0.181, -0.071, -0.038 };
   double X[] = { -0.015, 0.132 };
   int incX = -1;
   double Y[] = { -0.449, -0.219 };
   int incY = -1;
   double y_expected[] = { -0.036098, 9.27e-04 };
   cblas_dspmv(order, uplo, N, alpha, A, X, incX, beta, Y, incY);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(Y[i], y_expected[i], dbleps, "dspmv(case 1149)");
     }
   };
  };

}

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc


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

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

void freeDoubleValue(CGDouble* v){
}

void freeVectorValue(CGVector* data){
	free(data->data);
}

void freeMatrixValue(CGMatrix* data){
	free(data->data);
}

void freeNode(CGraph* graph, CGNode* node){
	//printf("freeing node %d\n", node->type);
	if(node->result != NULL){
		freeResultNode(node->result);
		free(node->result);
		node->result = NULL;
	}

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

	switch(node->type){
		case CGNT_CONSTANT:{
			switch(node->constant->type){
				case CGVT_DOUBLE:
				{
					free(node->constant->value);
					free(node->constant);

					break;
				}

				case CGVT_VECTOR:{
					CGVector* vec = (CGVector*)node->constant->value;
					freeVectorValue(vec);
					free(node->constant->value);
					free(node->constant);

					break;
				}

				case CGVT_MATRIX:{
					CGMatrix* mat = (CGMatrix*)node->constant->value;
					freeMatrixValue(mat);
					free(node->constant->value);
					free(node->constant);
					break;
				}
			}
			break;
		}
		case CGNT_VARIABLE:{
			free(node->var);
			break;
		}
		case CGNT_BINARY_OPERATION:
			free(node->bop);
			break;
		case CGNT_UNARY_OPERATION:
			free(node->uop);
			break;
		case CGNT_AXIS_BOUND_OPERATION:
			free(node->axop);
			break;
		case CGNT_GRAPH:
			free(node->graph);
			break;
		case CGNT_CROSS_ENTROPY_LOSS_FUNC:
			free(node->crossEntropyLoss);
			break;
	}

}

void freeResultNode(CGResultNode* node){
	if(node->error){
		free(node->error);
		free(node);
		return;
	}

	switch(node->type){
		case CGVT_DOUBLE:
			freeDoubleValue(node->value);
			break;
		case CGVT_VECTOR:
			freeVectorValue(node->value);
			break;
		case CGVT_MATRIX:
			freeMatrixValue(node->value);
			break;
	}

	free(node->value);
}

void freeGraph(CGraph* graph){
	if(graph == NULL)
		return;

	const char *key;

	map_iter_t iter = map_iter(&graph->vars);

	while ((key = map_next(&graph->vars, &iter))) {
		CGNode* node = *map_get(&graph->vars, key);
		if(node != NULL){
			freeNode(graph, node);
			free(node);
		}
	}

	map_deinit(&graph->vars);

	int i = 0;
	CGNode* node;

	vec_foreach(&graph->nodes, node, i) {
		freeNode(graph, node);
		free(node);
	}

	vec_deinit(&graph->nodes);
	// graph pointer must be freed elsewhere. in lua API we create a copy so we cannot free the parameter of this function.
}

/*
 * Lua lightuser data are automatically freed within the lua VM,
 * this version is lua friendly.
 */
void freeGraph_lua(CGraph* graph){
	if(graph == NULL)
		return;

	const char *key;

	map_deinit(&graph->vars);

	int i = 0;
	CGNode* node;

	vec_deinit(&graph->nodes);
	// graph pointer must be freed elsewhere. in lua API we create a copy so we cannot free the parameter of this function.
}

void freeNode_lua(CGraph* graph, CGNode* node){
	//printf("freeing node %d\n", node->type);

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

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

	switch(node->type){
		case CGNT_CONSTANT:{
			switch(node->constant->type){
				case CGVT_DOUBLE:
				{
					free(node->constant->value);
					free(node->constant);

					break;
				}

				case CGVT_VECTOR:{
					CGVector* vec = (CGVector*)node->constant->value;
					freeVectorValue(vec);
					free(node->constant->value);
					free(node->constant);

					break;
				}

				case CGVT_MATRIX:{
					CGMatrix* mat = (CGMatrix*)node->constant->value;
					freeMatrixValue(mat);
					free(node->constant->value);
					free(node->constant);
					break;
				}
			}
			break;
		}
		case CGNT_VARIABLE:{
			free(node->var);
			break;
		}
		case CGNT_BINARY_OPERATION:
			free(node->bop);
			break;
		case CGNT_UNARY_OPERATION:
			free(node->uop);
			break;
		case CGNT_AXIS_BOUND_OPERATION:
			free(node->axop);
			break;
		case CGNT_GRAPH:
			free(node->graph);
			break;
		case CGNT_CROSS_ENTROPY_LOSS_FUNC:
			free(node->crossEntropyLoss);
			break;
	}

}

max_stars_repo_name: Enehcruon/cgraph

max_issues_repo_name: Enehcruon/cgraph

max_forks_repo_name: praisethemoon/ccgraph


text:
#pragma once

#include "PropertyType.h"

#include <CesiumUtility/SpanHelper.h>

#include <gsl/span>

#include <cassert>
#include <cstddef>

namespace CesiumGltf {
template <typename ElementType> class MetadataArrayView {
public:
  MetadataArrayView() : _valueBuffer{} {}

  MetadataArrayView(const gsl::span<const std::byte>& buffer) noexcept
      : _valueBuffer{
            CesiumUtility::reintepretCastSpan<const ElementType>(buffer)} {}

  const ElementType& operator[](int64_t index) const noexcept {
    return _valueBuffer[index];
  }

  int64_t size() const noexcept {
    return static_cast<int64_t>(_valueBuffer.size());
  }

private:
  gsl::span<const ElementType> _valueBuffer;
};

template <> class MetadataArrayView<bool> {
public:
  MetadataArrayView() : _valueBuffer{}, _bitOffset{0}, _instanceCount{0} {}

  MetadataArrayView(
      const gsl::span<const std::byte>& buffer,
      int64_t bitOffset,
      int64_t instanceCount) noexcept
      : _valueBuffer{buffer},
        _bitOffset{bitOffset},
        _instanceCount{instanceCount} {}

  bool operator[](int64_t index) const noexcept {
    index += _bitOffset;
    const int64_t byteIndex = index / 8;
    const int64_t bitIndex = index % 8;
    const int bitValue =
        static_cast<int>(_valueBuffer[byteIndex] >> bitIndex) & 1;
    return bitValue == 1;
  }

  int64_t size() const noexcept { return _instanceCount; }

private:
  gsl::span<const std::byte> _valueBuffer;
  int64_t _bitOffset;
  int64_t _instanceCount;
};

template <> class MetadataArrayView<std::string_view> {
public:
  MetadataArrayView()
      : _valueBuffer{}, _offsetBuffer{}, _offsetType{}, _size{0} {}

  MetadataArrayView(
      const gsl::span<const std::byte>& buffer,
      const gsl::span<const std::byte>& offsetBuffer,
      PropertyType offsetType,
      int64_t size) noexcept
      : _valueBuffer{buffer},
        _offsetBuffer{offsetBuffer},
        _offsetType{offsetType},
        _size{size} {}

  std::string_view operator[](int64_t index) const noexcept {
    const size_t currentOffset =
        getOffsetFromOffsetBuffer(index, _offsetBuffer, _offsetType);
    const size_t nextOffset =
        getOffsetFromOffsetBuffer(index + 1, _offsetBuffer, _offsetType);
    return std::string_view(
        reinterpret_cast<const char*>(_valueBuffer.data() + currentOffset),
        (nextOffset - currentOffset));
  }

  int64_t size() const noexcept { return _size; }

private:
  static size_t getOffsetFromOffsetBuffer(
      size_t instance,
      const gsl::span<const std::byte>& offsetBuffer,
      PropertyType offsetType) noexcept {
    switch (offsetType) {
    case PropertyType::Uint8: {
      assert(instance < offsetBuffer.size() / sizeof(uint8_t));
      const uint8_t offset = *reinterpret_cast<const uint8_t*>(
          offsetBuffer.data() + instance * sizeof(uint8_t));
      return static_cast<size_t>(offset);
    }
    case PropertyType::Uint16: {
      assert(instance < offsetBuffer.size() / sizeof(uint16_t));
      const uint16_t offset = *reinterpret_cast<const uint16_t*>(
          offsetBuffer.data() + instance * sizeof(uint16_t));
      return static_cast<size_t>(offset);
    }
    case PropertyType::Uint32: {
      assert(instance < offsetBuffer.size() / sizeof(uint32_t));
      const uint32_t offset = *reinterpret_cast<const uint32_t*>(
          offsetBuffer.data() + instance * sizeof(uint32_t));
      return static_cast<size_t>(offset);
    }
    case PropertyType::Uint64: {
      assert(instance < offsetBuffer.size() / sizeof(uint64_t));
      const uint64_t offset = *reinterpret_cast<const uint64_t*>(
          offsetBuffer.data() + instance * sizeof(uint64_t));
      return static_cast<size_t>(offset);
    }
    default:
      assert(false && "Offset type has unknown type");
      return 0;
    }
  }
  gsl::span<const std::byte> _valueBuffer;
  gsl::span<const std::byte> _offsetBuffer;
  PropertyType _offsetType;
  int64_t _size;
};
} // namespace CesiumGltf

max_stars_repo_name: JiangMuWen/cesium-native

max_issues_repo_name: JiangMuWen/cesium-native

max_forks_repo_name: JiangMuWen/cesium-native


text:
/*
 *  iaf_chxk_2008.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 IAF_CHXK_2008_H
#define IAF_CHXK_2008_H

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

#ifdef HAVE_GSL

// C 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 "recordables_map.h"
#include "ring_buffer.h"
#include "universal_data_logger.h"

/* BeginDocumentation
Name: iaf_chxk_2008 - Conductance based leaky integrate-and-fire neuron model
                      used in Casti et al 2008.

Description:
iaf_chxk_2008 is an implementation of a spiking neuron using IAF dynamics with
conductance-based synapses [1]. It is modeled after iaf_cond_alpha with the
addition of after hyper-polarization current instead of a membrane potential
reset. Incoming spike events induce a post-synaptic change of conductance
modeled by an alpha function. The alpha function is normalized such that an
event of weight 1.0 results in a peak current of 1 nS at t = tau_syn.

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

V_m        double - Membrane potential in mV
E_L        double - Leak reversal potential in mV.
C_m        double - Capacity of the membrane in pF
V_th       double - Spike threshold in mV.
E_ex       double - Excitatory reversal potential in mV.
E_in       double - Inhibitory reversal potential in mV.
g_L        double - Leak conductance in nS.
tau_ex     double - Rise time of the excitatory synaptic alpha function in ms.
tau_in     double - Rise time of the inhibitory synaptic alpha function in ms.
I_e        double - Constant input current in pA.
tau_ahp    double - Afterhyperpolarization (AHP) time constant in ms.
E_ahp      double - AHP potential in mV.
g_ahp      double - AHP conductance in nS.
ahp_bug    bool   - Defaults to false. If true, behaves like original
                    model implementation.

References:
[1] Casti A, Hayot F, Xiao Y, and Kaplan E (2008) A simple model of retina-LGN
transmission. J Comput Neurosci 24:235-252.

Sends: SpikeEvent

Receives: SpikeEvent, CurrentEvent

Author: Heiberg

SeeAlso: iaf_cond_alpha
*/

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 iaf_chxk_2008_dynamics( double, const double*, double*, void* );

/**
 * Integrate-and-fire neuron model with two conductance-based synapses.
 */
class iaf_chxk_2008 : public Archiving_Node
{

  // Boilerplate function declarations --------------------------------

public:
  iaf_chxk_2008();
  iaf_chxk_2008( const iaf_chxk_2008& );
  ~iaf_chxk_2008();

  /**
   * 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 );

  bool
  is_off_grid() const
  {
    return true;
  } // uses off_grid events

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

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

  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 );

  // END Boilerplate function declarations ----------------------------

  // Friends --------------------------------------------------------

  // make dynamics function quasi-member
  friend int iaf_chxk_2008_dynamics( double, const double*, double*, void* );

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

private:
  // Parameters class -------------------------------------------------

  //! Model parameters
  struct Parameters_
  {
    double V_th; //!< Threshold Potential in mV
    double g_L;  //!< Leak Conductance in nS
    double C_m;  //!< Membrane Capacitance in pF
    double E_ex; //!< Excitatory reversal Potential in mV
    double E_in; //!< Inhibitory reversal Potential in mV
    double E_L;  //!< Leak reversal Potential (a.k.a. resting potential) in mV
    double tau_synE; //!< Synaptic Time Constant Excitatory Synapse in ms
    double tau_synI; //!< Synaptic Time Constant for Inhibitory Synapse in ms
    double I_e;      //!< Constant Current in pA
    double tau_ahp;  //!< Afterhyperpolarization (AHP) time constant
    double g_ahp;    //!< AHP conductance
    double E_ahp;    //!< AHP potential
    bool ahp_bug;    //!< If true, discard AHP conductance value from previous
                     //!< spikes
    Parameters_();   //!< Set default parameter values

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

  // State variables class --------------------------------------------

  /**
   * State variables of the model.
   *
   * State variables consist of the state vector for the subthreshold
   * dynamics and the refractory count. The state vector must be a
   * C-style array to be compatible with GSL ODE solvers.
   *
   * @note Copy constructor and assignment operator are required because
   *       of the C-style array.
   */
public:
  struct State_
  {

    //! Symbolic indices to the elements of the state vector y
    enum StateVecElems
    {
      V_M = 0,
      DG_EXC,
      G_EXC,
      DG_INH,
      G_INH,
      DG_AHP,
      G_AHP,
      STATE_VEC_SIZE
    };

    //! state vector, must be C-array for GSL solver
    double y[ STATE_VEC_SIZE ];

    //!< number of refractory steps remaining
    int r;

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

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

    /**
     * Set state from values in dictionary.
     * Requires Parameters_ as argument to, e.g., check bounds.'
     */
    void set( const DictionaryDatum&, const Parameters_& );
  };

private:
  // Buffers class --------------------------------------------------------

  /**
   * Buffers of the model.
   * Buffers are on par with state variables in terms of persistence,
   * i.e., initialized only upon first Simulate call after ResetKernel
   * or ResetNetwork, but are implementation details hidden from the user.
   */
  struct Buffers_
  {
    Buffers_( iaf_chxk_2008& );                  //!<Sets buffer pointers to 0
    Buffers_( const Buffers_&, iaf_chxk_2008& ); //!<Sets buffer pointers to 0

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

    /** buffers and sums up incoming spikes/currents */
    RingBuffer spike_exc_;
    RingBuffer spike_inh_;
    RingBuffer currents_;

    /* GSL ODE stuff */
    gsl_odeiv_step* s_;    //!< stepping function
    gsl_odeiv_control* c_; //!< adaptive step size 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_;            //!< step size in ms
    double IntegrationStep_; //!< current integration time step, updated by GSL

    /**
     * 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_;
  };

  // Variables class -------------------------------------------------------

  /**
   * Internal variables of the model.
   * Variables are re-initialized upon each call to Simulate.
   */
  struct Variables_
  {
    /**
     * Impulse to add to DG_EXC on spike arrival to evoke unit-amplitude
     * conductance excursion.
     */
    double PSConInit_E;

    /**
     * Impulse to add to DG_INH on spike arrival to evoke unit-amplitude
     * conductance excursion.
     */
    double PSConInit_I;

    /**
     * Impulse to add to DG_AHP on spike generation to evoke unit-amplitude
     * conductance excursion.
     */
    double PSConInit_AHP;
  };

  // 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 ];
  }

  //! Read out remaining refractory time, used by UniversalDataLogger
  double
  get_r_() const
  {
    return Time::get_resolution().get_ms() * S_.r;
  }

  double
  get_I_syn_exc_() const
  {
    return S_.y[ State_::G_EXC ] * ( S_.y[ State_::V_M ] - P_.E_ex );
  }
  double
  get_I_syn_inh_() const
  {
    return S_.y[ State_::G_INH ] * ( S_.y[ State_::V_M ] - P_.E_in );
  }
  double
  get_I_ahp_() const
  {
    return S_.y[ State_::G_AHP ] * ( S_.y[ State_::V_M ] - P_.E_ahp );
  }

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

  // keep the order of these lines, seems to give best performance
  Parameters_ P_;
  State_ S_;
  Variables_ V_;
  Buffers_ B_;

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

// Boilerplate inline function definitions ----------------------------------

inline port
iaf_chxk_2008::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
iaf_chxk_2008::handles_test_event( SpikeEvent&, rport receptor_type )
{
  if ( receptor_type != 0 )
  {
    throw UnknownReceptorType( receptor_type, get_name() );
  }
  return 0;
}

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

inline port
iaf_chxk_2008::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
iaf_chxk_2008::get_status( DictionaryDatum& d ) const
{
  P_.get( d );
  S_.get( d );
  Archiving_Node::get_status( d );

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

inline void
iaf_chxk_2008::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, ptmp );   // 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 // IAF_CHXK_2008_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:
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_math.h>
#include "gsl_cblas.h"

#include "tests.h"

void
test_syr2 (void) {
const double flteps = 1e-4, dbleps = 1e-6;
  {
   int order = 101;
   int uplo = 121;
   int N = 1;
   int lda = 1;
   float alpha = 0.0f;
   float A[] = { 0.862f };
   float X[] = { 0.823f };
   int incX = -1;
   float Y[] = { 0.699f };
   int incY = -1;
   float A_expected[] = { 0.862f };
   cblas_ssyr2(order, uplo, N, alpha, X, incX, Y, incY, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[i], A_expected[i], flteps, "ssyr2(case 1434)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int N = 1;
   int lda = 1;
   float alpha = 0.0f;
   float A[] = { 0.862f };
   float X[] = { 0.823f };
   int incX = -1;
   float Y[] = { 0.699f };
   int incY = -1;
   float A_expected[] = { 0.862f };
   cblas_ssyr2(order, uplo, N, alpha, X, incX, Y, incY, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[i], A_expected[i], flteps, "ssyr2(case 1435)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int N = 1;
   int lda = 1;
   float alpha = 0.0f;
   float A[] = { 0.862f };
   float X[] = { 0.823f };
   int incX = -1;
   float Y[] = { 0.699f };
   int incY = -1;
   float A_expected[] = { 0.862f };
   cblas_ssyr2(order, uplo, N, alpha, X, incX, Y, incY, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[i], A_expected[i], flteps, "ssyr2(case 1436)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int N = 1;
   int lda = 1;
   float alpha = 0.0f;
   float A[] = { 0.862f };
   float X[] = { 0.823f };
   int incX = -1;
   float Y[] = { 0.699f };
   int incY = -1;
   float A_expected[] = { 0.862f };
   cblas_ssyr2(order, uplo, N, alpha, X, incX, Y, incY, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[i], A_expected[i], flteps, "ssyr2(case 1437)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 121;
   int N = 1;
   int lda = 1;
   double alpha = 0;
   double A[] = { -0.824 };
   double X[] = { 0.684 };
   int incX = -1;
   double Y[] = { 0.965 };
   int incY = -1;
   double A_expected[] = { -0.824 };
   cblas_dsyr2(order, uplo, N, alpha, X, incX, Y, incY, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[i], A_expected[i], dbleps, "dsyr2(case 1438)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int N = 1;
   int lda = 1;
   double alpha = 0;
   double A[] = { -0.824 };
   double X[] = { 0.684 };
   int incX = -1;
   double Y[] = { 0.965 };
   int incY = -1;
   double A_expected[] = { -0.824 };
   cblas_dsyr2(order, uplo, N, alpha, X, incX, Y, incY, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[i], A_expected[i], dbleps, "dsyr2(case 1439)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int N = 1;
   int lda = 1;
   double alpha = 0;
   double A[] = { -0.824 };
   double X[] = { 0.684 };
   int incX = -1;
   double Y[] = { 0.965 };
   int incY = -1;
   double A_expected[] = { -0.824 };
   cblas_dsyr2(order, uplo, N, alpha, X, incX, Y, incY, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[i], A_expected[i], dbleps, "dsyr2(case 1440)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int N = 1;
   int lda = 1;
   double alpha = 0;
   double A[] = { -0.824 };
   double X[] = { 0.684 };
   int incX = -1;
   double Y[] = { 0.965 };
   int incY = -1;
   double A_expected[] = { -0.824 };
   cblas_dsyr2(order, uplo, N, alpha, X, incX, Y, incY, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[i], A_expected[i], dbleps, "dsyr2(case 1441)");
     }
   };
  };

}


text:
#include <math.h>
#include <string.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_sf_erf.h>
#include "prob_functions.h"
#include "data.h"

static const double ALPHA_PRIOR_MEAN = 1, BETA_PRIOR_MEAN = 1;

double my_f (const gsl_vector *x, void *params)
{
	Dataset *data = (Dataset *) params;

	unpackX(x, data);

	return - computeQ(data);
}

void my_df (const gsl_vector *x, void *params, gsl_vector *g)
{
	int i, j;

	Dataset *data = (Dataset *) params;
	double *dQdAlpha = (double *) malloc(sizeof(double) * data->numLabelers * data->numWorkflows);
	double *dQdBeta = (double *) malloc(sizeof(double) * data->numImages * data->numWorkflows);

	unpackX(x, data);

	gradientQ(data, dQdAlpha, dQdBeta);

	/* Pack dQdAlpha and dQdBeta into gsl_vector */
	for (i = 0; i < data->numLabelers * data->numWorkflows; i++) {
		gsl_vector_set(g, i, - dQdAlpha[i]);  /* Flip the sign since we want to minimize */
	}
	for (j = 0; j < data->numImages * data->numWorkflows; j++) {
	  gsl_vector_set(g, (data->numLabelers * data->numWorkflows) + j, - dQdBeta[j]);  /* Flip the sign since we want to minimize */
	}

	free(dQdAlpha);
	free(dQdBeta);
}

void my_fdf (const gsl_vector *x, void *params, double *f, gsl_vector *g)
{
	*f = my_f(x, params);
	my_df(x, params, g);
}

void gradientQ (Dataset *data, double *dQdAlpha, double *dQdBeta)
{
	int i, j;
	int idx;

	/* This comes from the priors */
	/*RandoPriors*/
	/*
	for (i = 0; i < data->numLabelers; i++) {
		dQdAlpha[i] = - 1 * (data->alpha[i] - data->priorAlpha[i]);
	}

	for (j = 0; j < data->numImages; j++) {
		dQdBeta[j] = - 5 * (data->beta[j] - data->priorBeta[j]);
	}
	*/
	/*Here's where the real shit starts*/
	for (idx = 0; idx < data->numLabels; idx++) {
		int i = data->labels[idx].labelerId;
		int j = data->labels[idx].imageIdx;
		int k = data->labels[idx].workflowId;
		int alphaidx = k * data->numLabelers + i;
		int betaidx = k * data->numImages + j;
		int lij = data->labels[idx].label;
		double sigma = prob(data->alpha[alphaidx], data->beta[betaidx]);

		double increa, increb;
		double weight = data->probZ1[j] * (lij - sigma) + data->probZ0[j] * (1 - lij - sigma);
		//double weight = data->probZ1[j] * (2 * lij - 1) + data->probZ0[j] * (-2 * lij + 1);
		increa = .5 * weight * log(1-data->beta[betaidx]) * pow( 1-data->beta[betaidx],
									 data->alpha[alphaidx] );
		increb = -.5 * weight * data->alpha[alphaidx] * pow( 1-data->beta[betaidx],
								     data->alpha[alphaidx]-1 );
		dQdAlpha[alphaidx] += increa;
		dQdBeta[betaidx] += increb;
		//printf("Alpha: %f Beta: %f", data->alpha[alphaidx], data->beta[betaidx]);
		//printf("Blah%f %f %f\n", log(1-data->beta[j]), 1-data->beta[j], data->alpha[i]);
		//printf("%d %d %d %f %f %f %f\n", i, j, lij, sigma, increa, increb, weight);
		//printf("%f %f\n", data->probZ1[j], data->probZ0[j]);
		/*
		if ( isnan(prob(data->alpha[i], data->beta[j])) )
{
printf("%d %d\n", i, j);
printf("%f %f\n", data->probZ1[j], data->probZ0[j]);
printf("alpha %f beta %f\n", data->alpha[i], data->beta[j]);
printf("sigma %f\n", sigma);
printf("label %d\n", lij);
printf("prob: %f\n", prob);
printf("%f\n", pow( 1-data->beta[j], data->alpha[i]-1 ));
printf("%f\n", data->alpha[i]);
printf("%f %f\n", increa, increb);
}
fflush(0);

	for (i = 0; i < data->numLabelers; i++) {
		printf("i %d %f %f\n", i, data->alpha[i], dQdAlpha[i]);
	}

	for (j = 0; j < data->numImages; j++) {
		printf("j %d %f %f\n", j, data->beta[j], dQdBeta[j]);
	}
		*/
	  }

}
/*
double computeLikelihood (Dataset *data)
{
	int i, j;
	int idx;
	double L = 0;

	double *alpha = data->alpha, *beta = data->beta;

	for (j = 0; j < data->numImages; j++) {
		double P1 = data->priorZ1[j];
		double P0 = 1 - data->priorZ1[j];
		for (idx = 0; idx < data->numLabels; idx++) {
			if (data->labels[idx].imageIdx == j) {
				int i = data->labels[idx].labelerId;
				int lij = data->labels[idx].label;
		double sigma = prob(alpha[i], beta[j]);
				P1 *= pow(sigma, lij) * pow(1 - sigma, 1 - lij);
				P0 *= pow(sigma, 1 - lij) * pow(1 - sigma, lij);
			}
		}
		L += log(P1 + P0);
	}

	/// Add Gaussian (standard normal) prior for alpha
	for (i = 0; i < data->numLabelers; i++) {
		L += log(gsl_sf_erf_Z(alpha[i] - data->priorAlpha[i]));
	}

	/// Add Gaussian (standard normal) prior for beta
	for (j = 0; j < data->numImages; j++) {
		L += log(gsl_sf_erf_Z(beta[j] - data->priorBeta[j]));
	}

	return L;
}*/

double computeQ (Dataset *data)
{
	int i, j;
	int idx;
	double Q = 0;

	double *alpha = data->alpha, *beta = data->beta;

	/* Start with the expectation of the sum of priors over all images */
	for (j = 0; j < data->numImages; j++) {
		Q += data->probZ1[j] * log(data->priorZ1[j]);
		Q += data->probZ0[j] * log(1 - data->priorZ1[j]);
	}

	for (idx = 0; idx < data->numLabels; idx++) {
		int i = data->labels[idx].labelerId;
		int j = data->labels[idx].imageIdx;
		int k = data->labels[idx].workflowId;
		int alphaidx = k * data->numLabelers + i;
		int betaidx = k * data->numImages + j;
		int lij = data->labels[idx].label;
		/* Do some analytic manipulation first for numerical stability! */
		double logSigma = log( prob(alpha[alphaidx], beta[betaidx]) );
		double logOneMinusSigma = log( 1-prob(alpha[alphaidx], beta[betaidx]) );
		if (beta[betaidx] > 1 || beta[betaidx] < 0 || alpha[alphaidx] < 0) {
		  Q += -100000000000000000;
		} else {
		  Q += data->probZ1[j] * (lij * logSigma + (1 - lij) * logOneMinusSigma) +
		    data->probZ0[j] * ((1 - lij) * logSigma + lij * logOneMinusSigma);

		}
		if (isnan(Q)) {
		  printf("%d %d %d %f %f %f %f %f %f %f %f\n", i, j, lij, logSigma,
			 logOneMinusSigma,
			 prob(alpha[i], beta[j]),
			 data->probZ1[j],
			 data->probZ0[j],
			 alpha[alphaidx],
			 beta[betaidx],
			 Q);
		  abort();
		}
	}

	/* Add Gaussian (standard normal) prior for alpha */
	for (i = 0; i < data->numLabelers * data->numWorkflows; i++) {
		//Q += log(gsl_sf_erf_Z(alpha[i] - data->priorAlpha[i]));
	}

	/* Add Gaussian (standard normal) prior for beta */
	for (j = 0; j < data->numImages * data->numWorkflows; j++) {
		//Q += log(gsl_sf_erf_Z(beta[j] - data->priorBeta[j]));
	}

	// add penalty of beta not between 0 and 1
	for (j = 0; j < data->numImages * data->numWorkflows; j++) {
	  if ( data->beta[j] > 1 || data->beta[j] < 0 ) {
	    Q -= 100;
	  }
	}
	return Q;
}

double logProbL (int l, int z, double alphaI, double betaJ)
{
	double p;

	if (z == l) {
		p = log(prob(alphaI, betaJ));
	} else {
		p = log(1-prob(alphaI, betaJ));
	}
/*
  printf("p: %f\n", prob(alphaI, betaJ) );
  printf("1-p: %f\n", 1-prob(alphaI, betaJ) );
  printf("log: %f\n", p);
*/
	return p;
}

double prob(double alpha, double beta)
{
  return .5 + .5 * pow(1-beta, alpha);
}

void EStep (Dataset *data)
{
	int j;
	int idx;
/*
  printf("begin Estep\n");
  fflush(0);
*/
	for (j = 0; j < data->numImages; j++) {
		data->probZ1[j] = log(data->priorZ1[j]);
		data->probZ0[j] = log(1 - data->priorZ1[j]);
		//printf("%d %f %f\n", j, data->probZ1[j], data->probZ0[j]);
	}

	for (idx = 0; idx < data->numLabels; idx++) {
		int i = data->labels[idx].labelerId;
		int j = data->labels[idx].imageIdx;
		int k = data->labels[idx].workflowId;
		int alphaidx = k * data->numLabelers + i;
		int betaidx = k * data->numImages + j;
		int lij = data->labels[idx].label;
		//printf("JOE%d %d %d %f %f\n", k, alphaidx, betaidx, data->alpha[alphaidx], data->beta[betaidx]);
		data->probZ1[j] += logProbL(lij, 1, data->alpha[alphaidx],
					    data->beta[betaidx]);
		data->probZ0[j] += logProbL(lij, 0, data->alpha[alphaidx],
					    data->beta[betaidx]);
		//printf("%d %d %f %f\n", j, lij, data->probZ1[j], data->probZ0[j]);
	}

	/* Exponentiate and renormalize */
	for (j = 0; j < data->numImages; j++) {
	  //printf("Before: %f %f\n", data->probZ1[j], data->probZ0[j]);
	  /*Chris adding code here to do this better */
	  double normalizingDenominator = data->probZ1[j] +
	    log(1 + exp(data->probZ0[j] -
			data->probZ1[j]));
	  data->probZ1[j] = data->probZ1[j] - normalizingDenominator;
	  data->probZ0[j] = data->probZ0[j] - normalizingDenominator;
	  //printf("Hmm: %f %f %f\n", normalizingDenominator, data->probZ1[j], data->probZ0[j]);
	  data->probZ1[j] = exp(data->probZ1[j]);
	  data->probZ0[j] = exp(data->probZ0[j]);
	  /*
	    data->probZ1[j] = exp(data->probZ1[j]);
	    data->probZ0[j] = exp(data->probZ0[j]);
	    printf("%f %f\n", data->probZ1[j], data->probZ0[j]);
	    data->probZ1[j] = data->probZ1[j] / (data->probZ1[j] + data->probZ0[j]);
	    data->probZ0[j] = 1 - data->probZ1[j];
	  */
	  //printf("%d %f %f\n", j, data->probZ1[j], data->probZ0[j]);
	  if (isnan(data->probZ1[j])) {
	    printf("ABORTING?");
	    abort();
	  }
	}
/*
  printf("end Estep\n");
  fflush(0);
*/
}

void MStep (Dataset *data)
{
	double lastF;
	int iter, status;
	int i, j;

	int numLabelers = data->numLabelers;
	int numImages = data->numImages;
	int numWorkflows = data->numWorkflows;

	gsl_vector *x;
	gsl_multimin_fdfminimizer *s;
	const gsl_multimin_fdfminimizer_type *T;
	gsl_multimin_function_fdf my_func;

	x = gsl_vector_alloc((numLabelers * numWorkflows) + (numImages * numWorkflows));

	/* Pack alpha and beta into a gsl_vector */
	packX(x, data);

	/* Initialize objective function */
	my_func.n = (numLabelers * numWorkflows) + (numImages * numWorkflows);
	my_func.f = &my_f;
	my_func.df = &my_df;
	my_func.fdf = &my_fdf;
	my_func.params = data;

	/* Initialize minimizer */
	T = gsl_multimin_fdfminimizer_conjugate_pr;
	s = gsl_multimin_fdfminimizer_alloc(T, (numLabelers * numWorkflows) + (numImages * numWorkflows));

	gsl_multimin_fdfminimizer_set(s, &my_func, x, 0.01, 0.01);
	iter = 0;
	do {
		lastF = s->f;
		iter++;
		/*printf("iter=%d f=%f\n", iter, s->f);*/

		status = gsl_multimin_fdfminimizer_iterate(s);
		if (status) {
			break;
		}
		status = gsl_multimin_test_gradient(s->gradient, 1e-3);
//	} while (lastF - s->f > 0.01 && status == GSL_CONTINUE && iter < 25);
	} while (iter < 5);

	/* Unpack alpha and beta from gsl_vector */
	unpackX(s->x, data);

	gsl_multimin_fdfminimizer_free(s);
	gsl_vector_free(x);
}

max_stars_repo_name: ShreyaR/WorkerPoolSelection

max_issues_repo_name: ShreyaR/WorkerPoolSelection

max_forks_repo_name: ShreyaR/WorkerPoolSelection


text:
#include <stdio.h>
#include <stdlib.h>

// include files for optimized libraries
#if defined USE_ESSL
#include <essl.h>
#elif defined USE_MKL
#include <mkl_cblas.h>
#include <mkl_lapacke.h>
#elif defined USE_LAPACK
#include <cblas.h>
#include <lapacke.h>
#endif

// interface to f2c code

#include "f2c.h"
#include "mrrr.h"

static inline int PCA_ssytd2(char uplo, int n, real *a, int lda, real *d, real *e, real *tau)
{
    int info;
    ssytd2_(&uplo, &n, a, &lda, d, e, tau, &info, 1);
    return info;
}

static inline int PCA_sstemr(char jobz, char range, int n, real *d, real *e, real vl, real vu,
    int il, int iu, int *m, real *w, real *z, int ldz, int nzc, int *isuppz, int *tryrac,
    real *work, int lwork, int *iwork, int liwork)
{
    int info;
    pca_sstemr__(&jobz, &range, &n, d, e, &vl, &vu, &il, &iu, m, w, z, &ldz, &nzc, isuppz,
        tryrac, work, &lwork, iwork, &liwork, &info, 1, 1);
    return info;
}

static inline int PCA_sorm2l(char side, char trans, int m, int n, int k, real *a, int lda, real *tau, real *c, int ldc, real *work)
{
    int info;
    sorm2l_(&side, &trans, &m, &n, &k, a, &lda, tau, c, &ldc, work, &info, 1, 1);
    return info;
}

#ifndef NO_PULP
#include "utils.h"
#include "hwTrace.h"
#endif

#if 1

#define ALLOC(t, v, s) t v[s];
#define FREE(v)

#else

#define ALLOC(t, v, s) t *v = malloc(sizeof(t) * (s));
#define FREE(v) free(v);

#endif

// PCA main routine
// input is a column-major matrix with a row for each sample and a column for each variable
// output is a column-major matrix with a row for each sample and a column for each component

void PCA_mrrr(int samples, int variables, float *input, int components, float *output)
{
    int lwork = 18 * variables;
    int liwork = 10 * variables;
    ALLOC(real, A, variables * variables);
    ALLOC(real, T, samples * variables);
    ALLOC(real, d, variables);
    ALLOC(real, e, variables);
    ALLOC(real, tau, variables);
    ALLOC(int, isuppz, variables * 2);
    ALLOC(real, w, variables);
    ALLOC(real, Z, variables * variables);
    ALLOC(real, work, lwork);
    ALLOC(int, iwork, liwork);

    // pulp_trace_kernel_declare(0, "kernel 0");
    // pulp_trace_kernel_start(0, 1);

    // compute and subtract mean
    for (int j = 0; j < variables; j++) {
        real mean = 0.0;
        #pragma omp parallel for reduction(+:mean)
        for (int i = 0; i < samples; i++)
            mean += input[j * samples + i];
        mean /= samples;
        #pragma omp parallel for
        for (int i = 0; i < samples; i++)
            T[j * samples + i] = input[j * samples + i] - mean;
    }

    // compute A=T^T*T
    for (int j = 0; j < variables; j++)
        for (int i = 0; i <= j; i++) {
            real dot = 0;
            #pragma omp parallel for reduction(+:dot)
            for (int k = 0; k < samples; k++)
                dot += T[j * samples + k] * T[i * samples + k];
            A[i + j * variables] = dot;
        }

    // tridiagonalization
#if defined USE_MKL || USE_LAPACK
    int info = LAPACKE_ssytrd(LAPACK_COL_MAJOR, 'U', variables, A, variables, d, e, tau);
#else
    int info = PCA_ssytd2('U', variables, A, variables, d, e, tau);
#endif
    if (info != 0) {
        printf("Error in SSYTRD/SSYTD2: %i\n", info);
        abort();
    }

    // compute eigenvalues
    int il = variables - components + 1, iu = variables, m, tryrac = 1;
    real vl = 0.0, vu = 0.0;
    info = PCA_sstemr('V', 'I', variables, d, e, vl, vu, il, iu, &m, w, Z, variables, variables,
        isuppz, &tryrac, work, lwork, iwork, liwork);
    if (info != 0) {
        printf("Error in SSTEMR: %i\n", info);
        abort();
    }
    printf("%d: ", m);
    for (int i = 0; i < m; i++) printf("%d ", (int)w[i]); printf("\n");

    // compute eigenvectors
#if defined USE_MKL || USE_LAPACK
    info = LAPACKE_sormtr(LAPACK_COL_MAJOR, 'L', 'U', 'N', variables, m, A, variables, tau,
        Z, variables);
#else
    info = PCA_sorm2l('L', 'N', variables - 1, m, variables - 1, A + variables, variables, tau,
        Z, variables, work);
#endif
    if (info != 0) {
        printf("Error in SORMTR/SORM2L: %i\n", info);
        abort();
    }

#if defined USE_ESSL || USE_MKL || USE_LAPACK
    cblas_sgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, samples, components, variables,
        1.0, T, samples, Z, variables, 0.0, output, samples);
#else
    #pragma omp parallel for
    for (int i = 0; i < samples; i++)
        for (int j = 0; j < components; j++) {
            real t = 0;
            for (int k = 0; k < variables; k++)
                t += T[i + k * samples] * Z[j * variables + k];
            output[i + j * samples] = t;
        }
#endif

    FREE(T);
    FREE(A);
    FREE(e);
    FREE(d);
    FREE(tau);
    FREE(w);
    FREE(Z);
    FREE(isuppz);
    FREE(work);
    FREE(iwork);
    // pulp_trace_kernel_stop(0, 1);
}

max_stars_repo_name: stmach/micro-benchmarks

max_issues_repo_name: stmach/micro-benchmarks

max_forks_repo_name: stmach/micro-benchmarks