text:
/*---
  Flow*: A Verification Tool for Cyber-Physical Systems.
  Authors: Xin Chen, Sriram Sankaranarayanan, and Erika Abraham.
  Email: Xin Chen <chenxin415@gmail.com> if you have questions or comments.

  The code is released as is under the GNU General Public License (GPL).
---*/

#ifndef MATRIX_H_
#define MATRIX_H_

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

namespace flowstar
{

class Real;
class Interval;
class UnivariatePolynomial;
class Polynomial;
class RowVector;
class ColVector;
class rMatrix;
class iMatrix;
class upMatrix;
class mpMatrix;
class iMatrix2;
class HybridSystem;
class TaylorModelVec;

// The matrix class is implemented based on the data structure of gsl matrix.
class Matrix
{
private:
	gsl_matrix *data;
public:
	Matrix();
	Matrix(const int m, const int n);	// Create an m x n matrix, all of the entries are 0.
	Matrix(const int n);				// Create an n x n matrix, all of the entries are 0.
	Matrix(const Matrix & A);
	~Matrix();

	double get(const int i, const int j) const;				// Get the entry at position [i,j].
	void set(const double v, const int i, const int j);		// Set A[i,j] = v.

	int rows() const;
	int cols() const;

	void row(RowVector & result, const int i) const;			// Return the (i+1)-st row.

	void sortColumns();		// Sort the columns by size in descending order.
	int rank() const;

	void neg(Matrix & result) const;
	void neg_assign();

	void inverse(Matrix & result) const;
	void inverse_assign();

	void transpose(Matrix & result) const;
	void svd(Matrix & U) const;

	void QR(Matrix & D);
	void QRfactor(Matrix & Q);

	void output(FILE *fp) const;

	Matrix & operator += (const Matrix & A);
	Matrix & operator -= (const Matrix & A);
	Matrix & operator *= (const Matrix & A);

	Matrix operator + (const Matrix & A) const;
	Matrix operator - (const Matrix & A) const;
	Matrix operator * (const Matrix & A) const;

	Matrix & operator = (const Matrix & A);
};

class RowVector
{
private:
	Matrix vec;
public:
	RowVector();
	RowVector(const int n);
	RowVector(const RowVector & v);
	~RowVector();

	void set(const double v, const int pos);
	double get(const int pos) const;
	int size() const;

	void transpose(ColVector & result) const;

	void neg(RowVector & result) const;
	void neg_assign();

	void dump(FILE *fp) const;

	double innerProd(const RowVector & v) const;
	double EuclideanNorm() const;
	void normalize();

	bool operator == (const RowVector & v) const;

	RowVector & operator += (const RowVector & v);
	RowVector & operator -= (const RowVector & v);
	RowVector operator + (const RowVector & v) const;
	RowVector operator - (const RowVector & v) const;

	RowVector & operator = (const RowVector & v);

	friend class ColVector;
};

class ColVector
{
private:
	Matrix vec;
public:
	ColVector();
	ColVector(const int n);
	ColVector(const ColVector & v);
	~ColVector();

	void set(const double v, const int pos);
	double get(const int pos) const;
	int size() const;

	void transpose(RowVector & result) const;

	void neg(ColVector & result) const;
	void neg_assign();

	void mul(ColVector & result, const Matrix & m) const;
	void mul_assign(const Matrix & m);

	ColVector & operator += (const ColVector & v);
	ColVector & operator -= (const ColVector & v);
	ColVector operator + (const ColVector & v) const;
	ColVector operator - (const ColVector & v) const;

	ColVector & operator = (const ColVector & v);

	friend class RowVector;
};

class bMatrix
{
protected:
	bool *data;
	int size1;
	int size2;

public:
	bMatrix();
	bMatrix(const int m, const int n);
	bMatrix(const bMatrix & B);
	~bMatrix();

	int rows() const;
	int cols() const;

	void output(FILE *fp) const;

	bMatrix & operator += (const bMatrix & B);
	bool * operator [] (const int i);

	bMatrix & operator = (const bMatrix & B);
};

class rMatrix
{
protected:
	Real *data;
	int size1;
	int size2;

public:
	rMatrix();
	rMatrix(const int m, const int n);	// zero matrix
	rMatrix(const int n);				// identity matrix
	rMatrix(const rMatrix & rmatrix);
	~rMatrix();

	int rows() const;
	int cols() const;

	void abs(rMatrix & result) const;
	void transpose(rMatrix & result) const;

	void add_RNDD(rMatrix & result, const rMatrix & rmatrix) const;
	void add_assign_RNDD(const rMatrix & rmatrix);
	void add_RNDU(rMatrix & result, const rMatrix & rmatrix) const;
	void add_assign_RNDU(const rMatrix & rmatrix);

	void mul_RNDD(rMatrix & result, const rMatrix & rmatrix) const;
	void mul_RNDU(rMatrix & result, const rMatrix & rmatrix) const;

	void output(FILE *fp) const;

	rMatrix & operator += (const rMatrix & B);
	rMatrix & operator -= (const rMatrix & B);
	rMatrix & operator *= (const rMatrix & B);
	rMatrix operator + (const rMatrix & B) const;
	rMatrix operator - (const rMatrix & B) const;
	rMatrix operator * (const rMatrix & B) const;

	rMatrix & operator *= (const Real & r);
	rMatrix & operator *= (const double d);

	rMatrix & operator = (const rMatrix & rmatrix);		// should always be the same precision

	Real * operator [] (const int i);

	friend void to_iMatrix2(iMatrix2 & result, const rMatrix & lo, const rMatrix & up);

	friend class iMatrix;
	friend class iMatrix2;
	friend class upMatrix;
};

// matrix for intervals
class iMatrix
{
protected:
	Interval *data;
	int size1;
	int size2;

public:
	iMatrix();
	iMatrix(const int m, const int n);	// zero matrix
	iMatrix(const int n);				// identity matrix
	iMatrix(const iMatrix & A);
	iMatrix(const rMatrix & A);
	iMatrix(const iMatrix2 & A);
	iMatrix(const std::vector<Interval> & box);
	iMatrix(const std::string & matlab_format);
	~iMatrix();

	void clear();

	void mul(iMatrix & result, iMatrix & A);

	int rows() const;
	int cols() const;

	bool isZero() const;

	void pow(iMatrix & result, const int order) const;
	void pow_assign(const int order);

	double max_norm() const;
	void max_norm(Real & norm) const;

	void transpose(iMatrix & result) const;
	void times_pars(mpMatrix & result) const;
	void center();

	void linearTrans(std::vector<Polynomial> & result, const std::vector<Polynomial> & polyVec) const;

	void right_scale_assign(const std::vector<Interval> & scalars);

	void to_iMatrix2(iMatrix2 & A) const;

	void output(FILE *fp) const;

	iMatrix & operator += (const iMatrix & A);
	iMatrix & operator += (const Real & rad);
	iMatrix & operator -= (const iMatrix & A);
	iMatrix & operator *= (const iMatrix & A);
	iMatrix & operator *= (const Interval & I);
	iMatrix & operator *= (const double c);
	iMatrix & operator /= (const Interval & I);
	iMatrix & operator /= (const double c);

	iMatrix operator + (const iMatrix & A) const;
	iMatrix operator - (const iMatrix & A) const;
	iMatrix operator * (const iMatrix & A) const;
	iMatrix operator * (const Interval & I) const;
	iMatrix operator * (const double c) const;
	upMatrix operator * (const upMatrix & upm) const;

	mpMatrix operator * (const mpMatrix & mpm) const;
	TaylorModelVec operator * (const TaylorModelVec & tmv) const;

	Interval * operator [] (const int i);
	iMatrix & operator = (const iMatrix & A);
	iMatrix & operator = (const rMatrix & A);
	iMatrix & operator = (const iMatrix2 & A);
	iMatrix & operator = (const std::vector<Interval> & box);

	friend class upMatrix;
	friend class mpMatrix;
	friend class rMatrix;
	friend class iMatrix2;
	friend class Zonotope;
};

class iMatrix2
{
public:
	rMatrix center;
	rMatrix radius;

public:
	iMatrix2();
	iMatrix2(const int m, const int n);	// zero matrix
	iMatrix2(const int n);				// identity matrix
	iMatrix2(const iMatrix2 & A);
	iMatrix2(const iMatrix & A);
	iMatrix2(const std::vector<Interval> & box);
	~iMatrix2();

	int rows() const;
	int cols() const;

	void to_iMatrix(iMatrix & A) const;
	void transpose(iMatrix2 & result) const;
	void mag(Interval & I, const int i, const int j);
	void mag(Real & r, const int i, const int j);
	void add_assign(const Interval & I, const int i, const int j);

	iMatrix2 & operator += (const iMatrix2 & A);
	iMatrix2 & operator += (const iMatrix & A);
	iMatrix2 & operator += (const Real & rad);
	iMatrix2 & operator *= (const iMatrix2 & A);
	iMatrix2 & operator *= (const Interval & I);

	iMatrix2 operator + (const iMatrix2 & A) const;
	iMatrix2 operator * (const iMatrix2 & A) const;
	iMatrix operator * (const iMatrix & A) const;
	upMatrix operator * (const upMatrix & upm) const;
	TaylorModelVec operator * (const TaylorModelVec & tmv) const;

	iMatrix2 & operator = (const iMatrix2 & A);

	void output(FILE *fp) const;

	friend void to_iMatrix2(iMatrix2 & result, const rMatrix & lo, const rMatrix & up);

	friend class upMatrix;
	friend class mpMatrix;
	friend class rMatrix;
	friend class iMatrix;
	friend class Zonotope;
	friend class Polynomial;
	friend class HybridSystem;
};

// matrix for univariate polynomials
class upMatrix
{
protected:
	UnivariatePolynomial *data;
	int size1;
	int size2;

public:
	upMatrix();
	upMatrix(const int m, const int n);
	upMatrix(const int n);
	upMatrix(const iMatrix & A);
	upMatrix(const iMatrix2 & A);
	upMatrix(const upMatrix & upm);
	~upMatrix();

	int rows() const;
	int cols() const;
	int degree() const;

	bool isZero() const;

	void intEval(iMatrix & result, const std::vector<Interval> & val_exp_table) const;	// interval evaluation based on the monomial form
	void intEval(iMatrix & result, const Interval & val) const;							// interval evaluation based on the Horner form

	void intEval(iMatrix2 & result, const std::vector<Interval> & val_exp_table) const;
	void intEval(iMatrix2 & result, const Interval & val) const;

	void integral();
	void times_x(const int order);
	void transpose(upMatrix & result) const;

	void times_pars(mpMatrix & result) const;

	void ctrunc(iMatrix & rem, const int order, const std::vector<Interval> & val_exp_table);
	void ctrunc(iMatrix & rem, const int order, const Interval & val);

	void ctrunc(iMatrix & rem1, iMatrix & rem2, const int order, const std::vector<Interval> & val1_exp_table, const std::vector<Interval> & val2_exp_table);
	void ctrunc(iMatrix & rem1, iMatrix & rem2, const int order, const Interval & val1, const Interval & val2);

	void ctrunc(const int order, const std::vector<Interval> & val_exp_table);
	void ctrunc(const int order, const Interval & val);

	void nctrunc(const int order);

	void substitute(upMatrix & result, const std::vector<UnivariatePolynomial> & t_exp_table) const;
	void substitute(upMatrix & result, const UnivariatePolynomial & t) const;

	void output(FILE *fp) const;

	void decompose(upMatrix & positive, upMatrix & negative, iMatrix2 & im2_rem) const;

	upMatrix & operator += (const upMatrix & upm);
	upMatrix & operator += (const iMatrix & A);
	upMatrix & operator += (const Real & rad);
	upMatrix & operator -= (const upMatrix & upm);
	upMatrix & operator -= (const iMatrix & A);
	upMatrix & operator *= (const upMatrix & upm);
	upMatrix & operator *= (const iMatrix & A);
	upMatrix & operator *= (const Interval & I);

	upMatrix operator + (const upMatrix & upm) const;
	upMatrix operator + (const iMatrix & A) const;
	upMatrix operator - (const upMatrix & upm) const;
	upMatrix operator - (const iMatrix & A) const;
	upMatrix operator * (const upMatrix & upm) const;
	upMatrix operator * (const iMatrix & A) const;
	upMatrix operator * (const iMatrix2 & A) const;
	upMatrix operator * (const Interval & I) const;
	mpMatrix operator * (const mpMatrix & mpm) const;

	UnivariatePolynomial * operator [] (const int i);

	upMatrix & operator = (const upMatrix & upm);

	friend class iMatrix;
	friend class iMatrix2;
};

// matrix for multivariate polynomials
class mpMatrix
{
protected:
	Polynomial *data;
	int size1;
	int size2;

public:
	mpMatrix();
	mpMatrix(const int m, const int n);
	mpMatrix(const int n);
	mpMatrix(const mpMatrix & mpm);
	~mpMatrix();

	int rows() const;
	int cols() const;

	void intEval(mpMatrix & result, const std::vector<Interval> & val_exp_table) const;
	void intEval(iMatrix & result, const std::vector<Interval> & domain) const;

	void output(FILE *fp, const std::vector<std::string> & varNames) const;

	mpMatrix & operator += (const mpMatrix & mpm);

	mpMatrix operator + (const mpMatrix & mpm) const;

	Polynomial * operator [] (const int i);
	mpMatrix & operator = (const mpMatrix & mpm);

	friend class iMatrix;
	friend class upMatrix;
};

class MatrixParseSetting
{
public:
	std::string strExpression;
	iMatrix result;

public:
	MatrixParseSetting();
	MatrixParseSetting(const MatrixParseSetting & setting);
	~MatrixParseSetting();

	MatrixParseSetting & operator = (const MatrixParseSetting & setting);
};

void to_iMatrix2(iMatrix2 & result, const rMatrix & lo, const rMatrix & up);

extern MatrixParseSetting matrixParseSetting;

}

void parse_Matrix();

#endif /* MATRIX_H_ */

ext: h


text:
#ifndef basis_dcab07d2_8e6a_4150_b6fb_1eb09f015f76_h
#define basis_dcab07d2_8e6a_4150_b6fb_1eb09f015f76_h

#include <gslib\error.h>
//#include <gslib\string.h>
#include <gslib\std.h>
#include <rathen\config.h>
#include <rathen\buildin.h>

__rathen_begin__

class __gs_novtable type abstract
{
public:
    virtual ~type() {}
    virtual const gchar* get_name() const = 0;
    virtual int get_size() const = 0;
    virtual type* get_oprret(const oprinfo& opr) = 0;
    /*
     * Commutative attribute, if an operator was commutative, then a series of transformation could be done
     * in the expression optimization, for example,
     * it would always prefer to put the variables to the left side of the operator, and the constants to the
     * right side, so that we could generate the code try to put the constants within the instructions.
     */
    virtual bool is_comutative(const oprinfo& opr) const
    {
        if (opr.opr == _t("="))
            return false;
        return true;
    }
    /* Complex attribute, a complex type always has a construction/destruction operation. */
    virtual bool is_complex() const { return false; }
};

class __gs_novtable com_type abstract:
    public type
{
public:
    virtual bool is_complex() const { return true; }
};

class type_manager
{
public:
    typedef std::pair<string, type*> type_pair;
    typedef unordered_map<string, type*> type_map;
    typedef type_map::iterator iterator;
    typedef type_map::const_iterator const_iterator;

public:
    static type_manager* get_singleton_ptr()
    {
        static type_manager inst;
        return &inst;
    }
    ~type_manager()
    {
        std::for_each(_tpmap.begin(), _tpmap.end(), [](const type_pair& tp) { delete tp.second; });
        _tpmap.clear();
    }
    type* find_type(const gchar* name) { return find_type(string(name)); }
    type* find_type(const gchar* name, int len) { return find_type(string(name, len)); }
    type* find_type(const string& name)
    {
        auto i = _tpmap.find(name);
        return i == _tpmap.end() ? 0 : i->second;
    }

private:
    type_map    _tpmap;

public:
    template<class _cst>
    bool regtype()
    {
        type* t = gs_new(_cst);
        if(!_reginner(t)) { gs_del(type, t); return false; }
        return true;
    }

private:
    type_manager();
    bool _reginner(type* t)
    {
        assert(t);
        if(find_type(t->get_name()))
            return false;
        _tpmap.insert(std::make_pair(string(t->get_name()), t));
        return true;
    }
};

#define _type_manager type_manager::get_singleton_ptr()

class int_type:
    public type
{
public:
    virtual const gchar* get_name() const { return _t("int"); }
    virtual int get_size() const { return 4; }
    virtual type* get_oprret(const oprinfo& opr) { return this; }   // TODO
};

class bool_type:
    public type
{
public:
    virtual const gchar* get_name() const { return _t("bool"); }
    virtual int get_size() const { return 4; }
    virtual type* get_oprret(const oprinfo& opr) { return this; }   // TODO
};

class void_type:
    public type
{
public:
    virtual const gchar* get_name() const { return _t("void"); }
    virtual int get_size() const { return 0; }
    virtual type* get_oprret(const oprinfo& opr) { return 0; }
};

class string_type:
    public com_type
{
public:
    virtual const gchar* get_name() const { return _t("string"); }
    virtual int get_size() const { return (int)sizeof(gs::string); }
    virtual type* get_oprret(const oprinfo& opr) { return this; }   // TODO
};

class __gs_novtable object abstract
{
public:
    enum
    {
        tag_data,
        tag_const,
        tag_block,
        tag_node,
        tag_scope,
        //tag_
        //...
    };

protected:
    string      _name;

public:
    object() {}
    virtual ~object() {}
    virtual void set_name(const gchar* name) { _name = name; }
    virtual uint get_tag() const = 0;
    virtual bool is_holder() const = 0;
    virtual const string& get_name() const { return _name; }
};

struct unikey
{
public:
    typedef gchar* vckey;
    typedef const gchar* cckey;

protected:
    cckey       _key;

public:
    unikey(): _key(0) {}
    unikey(const gchar* k): _key(k) {}
    unikey(const string& k): _key(k.c_str()) {}
    unikey(const object* obj): _key(obj->get_name().c_str()) {}
    const gchar* get_key() const { return _key; }
};

struct indexing
{
public:
    struct hash
    {
        size_t operator()(const unikey& k) const
        { return string_hash(k.get_key()); }
    };
    struct equal
    {
        bool operator()(const unikey& k1, const unikey& k2) const
        { return string_hash(k1.get_key()) == string_hash(k2.get_key()); }
    };
    typedef unordered_map<unikey, object*, hash, equal> unimap;
    typedef unimap::iterator iterator;
    typedef unimap::const_iterator const_iterator;

protected:
    unimap      _pairs;

public:
    indexing() {}
    ~indexing() {}
    bool add_value(object* ptr)
    {
        assert(ptr);
        if(_pairs.insert(std::make_pair(unikey(ptr), ptr)).second) {
            set_error(_t("udt insert dp failed, maybe an error in name mangling."));
            return false;
        }
        return true;
    }
    object* find_value(const gchar* name)
    {
        assert(name);
        iterator i = _pairs.find(unikey(name));
        return i == _pairs.end() ? 0 :i->second;
    }
    const object* find_value(const gchar* name) const
    {
        assert(name);
        const_iterator i = _pairs.find(unikey(name));
        return i == _pairs.end() ? 0 : i->second;
    }
};

// class __gnvt pass_object:
//     virtual public object
// {
// public:
//     virtual bool is_holder() const { return false; }
// };
//
// class __gnvt hold_object:
//     virtual public object
// {
// public:
//     virtual bool is_holder() const { return true; }
//     virtual void add_indexing(object* obj) { _indexing.add_value(obj); }
//     virtual indexing& get_indexing() { return _indexing; }
//
// protected:
//     indexing    _indexing;
//
// public:
//     object* find_object(const gchar* name) { return _indexing.find_value(name); }
//     const object* find_object(const gchar* name) const { return _indexing.find_value(name); }
// };

// class block:
//     public object
// {
// protected:
//     indexing    _indexing;
//
// public:
//     virtual void add_indexing(object* obj) { _indexing.add_value(obj); }
//     virtual indexing& get_indexing() { return _indexing; }
//
// public:
//     object* find_object(const gchar* name) { return _indexing.find_value(name); }
//     const object* find_object(const gchar* name) const { return _indexing.find_value(name); }
// };
//
// class data:
//     public object
// {
// protected:
//     type*       _type;
//
// public:
//     virtual tag get_tag() const { return tag_data; }
//     virtual void set_type(type* t) { _type = t; }
//     virtual type* get_type() const { return _type; }
// };
//
// class reference:
//     public object
// {
// protected:
//     string      _origin;
//
// public:
//     virtual tag get_tag() const { return tag_const; }
//     virtual void set_origin(const gchar* str, int len) { _origin.assign(str, len); }
//     virtual const gchar* get_origin() const { return _origin.c_str(); }
//
// public:
//     void set_origin(const gchar* str) { set_origin(str, strtool::length(str)); }
// };

__rathen_end__

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include "2d_array.h"
#include "const.h"
#include "input.h"
#include "utilities.h"
#include "misc.h"

/******************************************************************************
MODULE:  greenband_test

PURPOSE:  Multitemporal cloud, cloud shadow, & snow masks (global version)

RETURN VALUE:
Type = int
ERROR error out due to memory allocation
SUCCESS no error encounted

HISTORY:
Date        Programmer       Reason
--------    ---------------  -------------------------------------
05282019    Su Ye

NOTES:
******************************************************************************/
int greenband_test
(
    int *clrx,
    float **clry,
    int start,
    int end,
    float rmse,
    float n_t,
    int *bl_ids,
    float *C0,
    float *C1
)
{
    char FUNC_NAME[] = "greenband_test";
    int i;
    float **x;
    float pred;
    int nums;
    float coefs[ROBUST_COEFFS];

    nums = end - start + 1;
    /* Allocate memory */
    x = (float **)allocate_2d_array(nums, ROBUST_COEFFS - 1, sizeof(float));
    if (x == NULL)
    {
        RETURN_ERROR("ERROR allocating x memory", FUNC_NAME, ERROR);
    }

    for (i = 0; i < nums; i++)
    {
        x[i][0] = (float)clrx[i+start];
    }

    /******************************************************************/
    /*                                                                */
    /* Do robust fitting for band 2 */
    /*                                                                */
    /******************************************************************/

    auto_robust_fit(x, clry, nums, start, 1, coefs);

    *C0 = coefs[0];
    *C1 = coefs[1];

    /******************************************************************/
    /*                                                                */
    /* predict band 2 and band 4 refs, bl_ids value of 0 is clear and   */
    /* 1 otherwise                                                    */
    /*                                                                */
    /******************************************************************/

    for (i = 0; i < nums; i++)
    {
        pred = coefs[0] + coefs[1] * (float)clrx[i+start];
        if (clry[1][i+start]-pred > (n_t * rmse))
        {
            // int testy = clry[1][i+start];
            // int testx = clrx[i+start];
            bl_ids[i] = 1;
        }
        else
        {
            bl_ids[i] = 0;
        }
    }

    /* Free allocated memory */
    if (free_2d_array ((void **) x) != SUCCESS)
    {
        RETURN_ERROR ("Freeing memory: x\n", FUNC_NAME, ERROR);
    }
    x = NULL;

    return (SUCCESS);
}

/******************************************************************************
MODULE:  nirband_test

PURPOSE:  nirband_test used to filter out shadow

RETURN VALUE:
Type = int
ERROR error out due to memory allocation
SUCCESS no error encounted

HISTORY:
Date        Programmer       Reason
--------    ---------------  -------------------------------------
05282019    Su Ye

NOTES:
******************************************************************************/
int nirband_test
(
    int *clrx,
    float **clry,
    int start,
    int end,
    float rmse,
    float n_t,
    int *bl_ids,
    float *C0,
    float *C1
)
{
    char FUNC_NAME[] = "nirband_test";
    int i;
    float **x;
    float pred;
    int nums;
    float coefs[ROBUST_COEFFS];

    nums = end - start + 1;
    /* Allocate memory */
    x = (float **)allocate_2d_array(nums, ROBUST_COEFFS - 1, sizeof(float));
    if (x == NULL)
    {
        RETURN_ERROR("ERROR allocating x memory", FUNC_NAME, ERROR);
    }

    for (i = 0; i < nums; i++)
    {
        x[i][0] = (float)clrx[i+start];
    }

    /******************************************************************/
    /*                                                                */
    /* Do robust fitting for band 4 */
    /*                                                                */
    /******************************************************************/

    auto_robust_fit(x, clry, nums, start, 3, coefs);

    *C0 = coefs[0];
    *C1 = coefs[1];

    /******************************************************************/
    /*                                                                */
    /* predict band 2 and band 4 refs, bl_ids value of 0 is clear and   */
    /* 1 otherwise                                                    */
    /*                                                                */
    /******************************************************************/

    for (i = 0; i < nums; i++)
    {
        pred = coefs[0] + coefs[1] * (float)clrx[i+start];
        if (clry[3][i+start]-pred < -(n_t * rmse))
        {
            bl_ids[i] = 1;
        }
        else
        {
            bl_ids[i] = 0;
        }
    }

    /* Free allocated memory */
    if (free_2d_array ((void **) x) != SUCCESS)
    {
        RETURN_ERROR ("Freeing memory: x\n", FUNC_NAME, ERROR);
    }
    x = NULL;

    return (SUCCESS);
}

/******************************************************************************
MODULE:  average_compositing

PURPOSE:  Running average compositing for pixel-based time series

RETURN VALUE:
Type = int (SUCCESS, ERROR or FAILURE)

HISTORY:
Date        Programmer       Reason
--------    ---------------  -------------------------------------
06/02/2019   Su Ye         Original Development
******************************************************************************/

int average_compositing
(
    short int **buf,            /* I:  pixel-based time series           */
    int *valid_date_array,      /* I: valid date time series               */
    int valid_date_count,       /* I: the number of valid dates               */
    int lower_ordinal,
    int upper_ordinal,
    int i_col,
    short int **out_compositing           /* O: outputted compositing results for four bands */
)
{
    int i, j;
    double index_sum[TOTAL_IMAGE_BANDS];
    int valid_count_window= 0;
    for(i = 0; i < TOTAL_IMAGE_BANDS; i++)
        index_sum[i] = 0;

    for(i = 0; i < valid_date_count; i++)
    {
        if((valid_date_array[i] > lower_ordinal - 1) && (valid_date_array[i] < upper_ordinal + 1))
        {
            for(j = 0; j < TOTAL_IMAGE_BANDS; j++)
            {
                index_sum[j] = index_sum[j] +  buf[j][i];
            }
            valid_count_window++;

        }

    }

    if(valid_count_window==0)
    {
        for(i = 0; i < TOTAL_IMAGE_BANDS; i++)
        {
            out_compositing[i][i_col] = -9999;
        }
        return SUCCESS;
    }

    for(j = 0; j < TOTAL_IMAGE_BANDS; j++)
    {
        out_compositing[j][i_col] = (short int)(index_sum[j] / valid_count_window);
    }
}

/******************************************************************************
MODULE:  average_compositing

PURPOSE:  Running average compositing for pixel-based time series

RETURN VALUE:
Type = int (SUCCESS, ERROR or FAILURE)

HISTORY:
Date        Programmer       Reason
--------    ---------------  -------------------------------------
06/03/2019   Su Ye         Original Development
******************************************************************************/

int median_compositing
(
    short int **buf,            /* I:  pixel-based time series           */
    int *valid_date_array,      /* I: valid date time series               */
    int valid_date_count,       /* I: the number of valid dates               */
    int i_col,
    short int **out_compositing           /* O: outputted compositing results for four bands */
)
{
    char FUNC_NAME[] = "median_compositing";
    short int *var;         /* pointer for allocation variable memory            */
    int i, j, m;

    if (valid_date_count == 1)
    {
        for (i = 0; i < TOTAL_IMAGE_BANDS; i++)
          out_compositing[i][i_col] = buf[i][0];
        return SUCCESS;
    }

    var = malloc(valid_date_count * sizeof(short int));
    if (var == NULL)
    {
        RETURN_ERROR ("Allocating var memory", FUNC_NAME, ERROR);
    }

    for (i = 0; i < TOTAL_IMAGE_BANDS; i++)
    {
        for (j = 0; j < valid_date_count; j++)
        {
            var[j] = buf[i][j];

        }
        quick_sort_float(var, 0, valid_date_count-1);
//            for (j = 0; j < dim2_end; j++)
//            {
//               printf("%f\n", var[j]);
//            }
        m = (valid_date_count) / 2;
        if (valid_date_count % 2 == 0)
        {
            //printf("%f\n", var[m-1]);
           //printf("%f\n", var[m]);
            out_compositing[i][i_col] = (short int)(var[m-1] + var[m]) / 2.0;
        }
        else
            out_compositing[i][i_col] = var[m];

    }

    free(var);
    return SUCCESS;

}

/******************************************************************************
MODULE:  hot_compositing

PURPOSE:  Running compositing for pixel-based time series

RETURN VALUE:
Type = int (SUCCESS, ERROR or FAILURE)

HISTORY:
Date        Programmer       Reason
--------    ---------------  -------------------------------------
05/02/2019   Su Ye         Original Development
******************************************************************************/
int hot_compositing
(
    short int **buf,            /* I:  pixel-based time series           */
    int *valid_date_array,      /* I: valid date time series               */
    int valid_date_count,       /* I: the number of valid dates               */
    int lower_ordinal,
    int upper_ordinal,
    int i_col,
    short int **out_compositing           /* O: outputted compositing results for four bands */
)
{
    int i, j;
    double wt;
    double index_sum[TOTAL_IMAGE_BANDS];
    for(i = 0; i < TOTAL_IMAGE_BANDS; i++)
        index_sum[i] = 0;
    double wt_sum = 0;
    int valid_count_window = 0;

    for(i = 0; i < valid_date_count; i++)
    {
        if((valid_date_array[i] > lower_ordinal - 1) && (valid_date_array[i] < upper_ordinal + 1))
        {
            wt = (double)1.0/((buf[BLUE_INDEX][i] - 0.5 * buf[RED_INDEX][i]) * (buf[BLUE_INDEX][i] - 0.5 * buf[RED_INDEX][i]));
            //wt = (float)1/(abs(buf[BLUE_INDEX][i] - 0.5 * buf[RED_INDEX][i]));
            for(j = 0; j < TOTAL_IMAGE_BANDS; j++)
            {
                index_sum[j] = index_sum[j] + buf[j][i] * wt;
            }
            wt_sum = wt_sum + wt;
            valid_count_window++;
        }
    }

    /********************************************/
    /*    condition 1: zero valid observation   */
    /********************************************/
    if(valid_count_window==0)
    {
        for(i = 0; i < TOTAL_IMAGE_BANDS; i++)
        {
            out_compositing[i][i_col] = -9999;
        }
        return SUCCESS;
    }

    /********************************************/
    /*    condition 2: standard procedures      */
    /********************************************/

    for(j = 0; j < TOTAL_IMAGE_BANDS; j++)
    {
        out_compositing[j][i_col] = (short int)(index_sum[j] / wt_sum);
        //out_compositing[j][i_col] = (short int)valid_count_window;
    }

    return SUCCESS;

}

/******************************************************************************
MODULE:  modified_hot_compositing

PURPOSE:  Running compositing for pixel-based time series by adding shadow consideration

RETURN VALUE:
Type = int (SUCCESS, ERROR or FAILURE)

HISTORY:
Date        Programmer       Reason
--------    ---------------  -------------------------------------
06/22/2019   Su Ye         Original Development
******************************************************************************/
int modified_hot_compositing
(
    short int **buf,            /* I:  pixel-based time series           */
    int *valid_date_array,      /* I: valid date time series               */
    int valid_date_count,       /* I: the number of valid dates               */
    int lower_ordinal,
    int upper_ordinal,
    int i_col,
    short int **out_compositing           /* O: outputted compositing results for four bands */
)
{
    int i, j;
    int status;
    double wt;
    double wt_shadow;
    double wt_cloud;
    double index_sum[TOTAL_IMAGE_BANDS];
    for(i = 0; i < TOTAL_IMAGE_BANDS; i++)
        index_sum[i] = 0;
    double wt_sum = 0;
    int valid_count_window = 0;
    short int** ts_subset;
    short int* ts_subset_selected_shadow;
    short int* ts_subset_selected_blue;
    short int variogram_shadow;
    short int medium_shadow;

    char FUNC_NAME[] = "modified_hot_compositing";

    ts_subset = (short int**)allocate_2d_array(TOTAL_IMAGE_BANDS, valid_date_count, sizeof(short int));
    if(ts_subset == NULL)
    {
        RETURN_ERROR ("Allocating ts_subset memory", FUNC_NAME, ERROR);
    }

    ts_subset_selected_shadow = (short int*)malloc(valid_date_count*sizeof(short int));
    if(ts_subset_selected_shadow == NULL)
    {
        RETURN_ERROR ("Allocating ts_subset_selected_shadow memory", FUNC_NAME, ERROR);
    }

    ts_subset_selected_blue = (short int*)malloc(valid_date_count*sizeof(short int));
    if(ts_subset_selected_blue == NULL)
    {
        RETURN_ERROR ("Allocating ts_subset_selected_blue memory", FUNC_NAME, ERROR);
    }

    for(i = 0; i < valid_date_count; i++)
    {
        if((valid_date_array[i] > lower_ordinal - 1) && (valid_date_array[i] < upper_ordinal + 1))
        {
            for(j = 0; j < TOTAL_IMAGE_BANDS; j++)
            {
                ts_subset[j][valid_count_window] = buf[j][i];
                if(j == NIR_INDEX)
                {
                   //ts_subset_selected_shadow[valid_count_window] = (buf[BLUE_INDEX][i] + buf[RED_INDEX][i] + buf[GREEN_INDEX][i])/3;
		   ts_subset_selected_shadow[valid_count_window] = buf[NIR_INDEX][i];
                    //printf("%i\n", ts_subset_selected[valid_count_window]);
                }
                if(j == BLUE_INDEX)
                {
                    ts_subset_selected_blue[valid_count_window] = buf[BLUE_INDEX][i];
                    //printf("%i\n", ts_subset_selected[valid_count_window]);
                }
            }

            valid_count_window++;

        }
    }

    /********************************************/
    /*    condition 1: zero valid observation   */
    /********************************************/
    if(valid_count_window < 3)
    {
        for(i = 0; i < TOTAL_IMAGE_BANDS; i++)
        {
            out_compositing[i][i_col] = -9999;
        }
        return SUCCESS;
    }

    //single_mean_rmse(ts_subset_selected, 0, valid_count_window - 1, &rmse, &mean);
    //quick_sort_shortint(ts_subset_selected, 0, valid_count_window - 1);

    //m = valid_count_window / 2;

    //medium_shadow = (ts_subset_selected[m] + ts_subset_selected[m - 1] + ts_subset_selected[m + 1])/3;
    single_median_variogram(ts_subset_selected_shadow, 0, valid_count_window - 1, &variogram_shadow, &medium_shadow);
    //single_median_quantile(ts_subset_selected_blue, 0, valid_count_window - 1, &quantile_blue, &medium_blue);
    //single_median_variogram(ts_subset_selected_hot, 0, valid_count_window - 1, &variogram_hot, &medium_hot);
    //penalty_slope = - 9.0 / (1000.0 * variogram);
//    penalty_slope = - 33.0 / (50.0 * variogram);
//    penalty_intercept = PENALTY_INTERCEPT;

//    for(i = 0; i < valid_count_window; i++)
//    {

//        if (ts_subset_selected[i] - medium < - 1.5 * variogram)
//        {
//            ts_subset_id[i] = 1;
//        }
//        else if (ts_subset_selected[i] - medium < - 2 * variogram)
//        {
//            ts_subset_id[i] = 2;
//        }
//        else if (ts_subset_selected[i] - medium < - 3 * variogram)
//        {
//            ts_subset_id[i] = 3;
//        }
//        else
//        {
//            ts_subset_id[i] = 0;
//        }

//    }
//    for(i = 0; i < valid_count_window; i++)
//    {

//            if (ts_subset_selected[i] - medium < - 1.5 * variogram)
//            {
//                ts_subset_id[i] = 1;
//            }
//            else
//            {
//                ts_subset_id[i] = 0;
//            }

//    }

    /********************************************/
    /*    condition 2: standard procedures      */
    /********************************************/
    double ratio;
    for(i = 0; i < valid_count_window; i++)
    {
         wt_cloud = (double) 1.0 / (ts_subset_selected_blue[i] * ts_subset_selected_blue[i]);

         if (ts_subset_selected_shadow[i] < medium_shadow)
         {
             ratio = (double)ts_subset_selected_shadow[i] / medium_shadow;
             wt_shadow = ratio * ratio * ratio * ratio;
         }
         else
             wt_shadow = 1.0;

         wt = wt_cloud * wt_shadow;
        //wt = (float)1/(abs(buf[BLUE_INDEX][i] - 0.5 * buf[RED_INDEX][i]));
        for(j = 0; j < TOTAL_IMAGE_BANDS; j++)
        {
            index_sum[j] = index_sum[j] + ts_subset[j][i] * wt;
        }
        wt_sum = wt_sum + wt;
    }

    for(j = 0; j < TOTAL_IMAGE_BANDS; j++)
    {
        out_compositing[j][i_col] = (short int)(index_sum[j] / wt_sum);
        //out_compositing[j][i_col] = (short int)valid_count_window;
    }

    /* free memory*/
    status = free_2d_array(ts_subset);
    if (status != SUCCESS)
    {
        RETURN_ERROR ("Freeing memory: ts_subset\n",
                      FUNC_NAME, FAILURE);
    }
    free(ts_subset_selected_shadow);
    free(ts_subset_selected_blue);

    return SUCCESS;

}

/******************************************************************************
MODULE:  valid_obs_count

PURPOSE:  count valid observation for each pixels

RETURN VALUE:
Type = int (SUCCESS, ERROR or FAILURE)

HISTORY:
Date        Programmer       Reason
--------    ---------------  -------------------------------------
06/22/2019   Su Ye         Original Development
******************************************************************************/
int valid_obs_count
(
    short int **buf,            /* I:  pixel-based time series           */
    int *valid_date_array,      /* I: valid date time series               */
    int valid_date_count,       /* I: the number of valid dates               */
    int lower_ordinal,
    int upper_ordinal,
    int i_col,
    short int **out_compositing           /* O: outputted compositing results for four bands */
)
{
    int index_sum[TOTAL_IMAGE_BANDS];
    int i, j, status;
    for(i = 0; i < TOTAL_IMAGE_BANDS; i++)
        index_sum[i] = 0;
    int valid_count_window = 0;
    short int** ts_subset;
    short int* ts_subset_selected;

    char FUNC_NAME[] = "modified_hot_compositing";

    ts_subset = (short int**)allocate_2d_array(TOTAL_IMAGE_BANDS, valid_date_count, sizeof(short int));
    if(ts_subset == NULL)
    {
        RETURN_ERROR ("Allocating ts_subset memory", FUNC_NAME, ERROR);
    }

    ts_subset_selected = (short int*)malloc(valid_date_count*sizeof(short int));
    if(ts_subset_selected == NULL)
    {
        RETURN_ERROR ("Allocating ts_subset_selected memory", FUNC_NAME, ERROR);
    }

    for(i = 0; i < valid_date_count; i++)
    {
        if((valid_date_array[i] > lower_ordinal - 1) && (valid_date_array[i] < upper_ordinal + 1))
        {
            for(j = 0; j < TOTAL_IMAGE_BANDS; j++)
            {
                ts_subset[j][valid_count_window] = buf[j][i];
                if(j == NIR_INDEX)
                {
                    ts_subset_selected[valid_count_window] = buf[NIR_INDEX][i];
                    //printf("%i\n", ts_subset_selected[valid_count_window]);
                }
            }
            valid_count_window++;

        }
    }

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

        out_compositing[j][i_col] = (short int)valid_count_window;
    }

    status = free_2d_array((void **)ts_subset);
    if (status != SUCCESS)
    {
        RETURN_ERROR ("Freeing memory: ts_subset\n",
                      FUNC_NAME, FAILURE);
    }
    free(ts_subset_selected);

    return SUCCESS;

}

/******************************************************************************
MODULE:  medium_compositing

PURPOSE:  Running compositing for pixel-based time series by adding shadow consideration

RETURN VALUE:
Type = int (SUCCESS, ERROR or FAILURE)

HISTORY:
Date        Programmer       Reason
--------    ---------------  -------------------------------------
06/22/2019   Su Ye         Original Development
******************************************************************************/
int medium_compositing
(
    short int **buf,            /* I:  pixel-based time series           */
    int *valid_date_array,      /* I: valid date time series               */
    int valid_date_count,       /* I: the number of valid dates               */
    int lower_ordinal,
    int upper_ordinal,
    int i_col,
    short int **out_compositing           /* O: outputted compositing results for four bands */
)
{
    int i, j, m;
    int status;
    double wt;
    double index_sum[TOTAL_IMAGE_BANDS];
    for(i = 0; i < TOTAL_IMAGE_BANDS; i++)
        index_sum[i] = 0;
    double wt_sum = 0;
    int valid_count_window = 0;
    short int** ts_subset;
    short int* ts_subset_selected;
    int* ts_subset_selected_index;
    char FUNC_NAME[] = "medium_compositing";

    ts_subset = (short int*)allocate_2d_array(TOTAL_IMAGE_BANDS, valid_date_count, sizeof(short int));
    if(ts_subset == NULL)
    {
        RETURN_ERROR ("Allocating ts_subset memory", FUNC_NAME, ERROR);
    }

    for(i = 0; i < valid_date_count; i++)
    {
        if((valid_date_array[i] > lower_ordinal - 1) && (valid_date_array[i] < upper_ordinal + 1))
        {
            for(j = 0; j < TOTAL_IMAGE_BANDS; j++)
            {
                ts_subset[j][valid_count_window] = buf[j][i];
            }
            valid_count_window++;
        }
    }

    ts_subset_selected = (short int*)malloc(valid_count_window*sizeof(short int));
    if(ts_subset_selected == NULL)
    {
        RETURN_ERROR ("Allocating ts_subset_selected memory", FUNC_NAME, ERROR);
    }

    ts_subset_selected_index = (int*)malloc(valid_count_window*sizeof(int));
    if(ts_subset_selected == NULL)
    {
        RETURN_ERROR ("Allocating ts_subset_selected_index memory", FUNC_NAME, ERROR);
    }

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

        ts_subset_selected[i] = ts_subset[NIR_INDEX][i];
        ts_subset_selected_index[i] = i;
            //printf("%i\n", ts_subset_selected[valid_count_window]);
    }

    /********************************************/
    /*    condition 1: zero valid observation   */
    /********************************************/
    if(valid_count_window == 0)
    {
        for(i = 0; i < TOTAL_IMAGE_BANDS; i++)
        {
            out_compositing[i][i_col] = -9999;
        }
        return SUCCESS;
    }

    quick_sort_shortint_index(ts_subset_selected, ts_subset_selected_index,
                              0, valid_count_window - 1);
//    for (j = 0; j < valid_count_window; j++)
//    {
//       printf("%i  %i\n", ts_subset_selected[j], ts_subset_selected_index[j]);
//    }

     m = valid_count_window / 2;
    if (valid_count_window % 2 == 0)
    {
        //printf("%f\n", var[m-1]);
       //printf("%f\n", var[m]);
        for(i = 0; i < TOTAL_IMAGE_BANDS; i++)
        {
           out_compositing[i][i_col] = (short int)((ts_subset[i][ts_subset_selected_index[m-1]] +
                   ts_subset[i][ts_subset_selected_index[m]]) / 2);
        }

    }
    else
    {
        for(i = 0; i < TOTAL_IMAGE_BANDS; i++)
        {
           out_compositing[i][i_col] = (short int)(ts_subset[i][ts_subset_selected_index[m]]);
        }
        //printf("%i\n", out_compositing[i][i_col]);
    }

    status = free_2d_array((void**)ts_subset);
    if (status != SUCCESS)
    {
        RETURN_ERROR ("Freeing memory: ts_subset\n",
                      FUNC_NAME, FAILURE);
    }
    free(ts_subset_selected);

    return SUCCESS;

}

/******************************************************************************
MODULE:  compositing_scanline

PURPOSE:  Running compositing for scanline-based time series

RETURN VALUE:
Type = int (SUCCESS, ERROR or FAILURE)

HISTORY:
Date        Programmer       Reason
--------    ---------------  -------------------------------------
05/02/2019   Su Ye         Original Development
******************************************************************************/
int compositing_scanline
(
    short int **buf,            /* I:  scanline-based time series           */
    int **valid_datearray_scanline,      /* I: valid date time series               */
    int *valid_datecount_scanline,       /* I: the number of valid dates               */
    int lower_ordinal,                /* I: lower ordinal date               */
    int upper_ordinal,                   /* I: upper_ordinal for temporal range of composition   */
    int num_samples,                /* I: the pixel number in a row              */
    int num_scenes,                 /* I: the number of scenes               */
    short int **out_compositing,           /* O: outputted compositing results for four bands */
    int method                       /* I: the compositing method{1 - fitting-weighted; 2 - fitting-normal; 3 - hot; 4 - average}  */
)
{
    int  j;
    int i_col;
    short int **tmp_buf;                   /* This is the image bands buffer, valid pixel only*/
    char FUNC_NAME[] = "compositing_scanline";
    int b_diagnosis = FALSE;
    Output_t* rec_c;

    tmp_buf = (short int **) allocate_2d_array (TOTAL_IMAGE_BANDS, num_scenes, sizeof (short int));
    if(tmp_buf == NULL)
    {
        RETURN_ERROR("ERROR allocating tmp_buf memory", FUNC_NAME, FAILURE);
    }

    rec_c = malloc(sizeof(Output_t));
    if(rec_c == NULL)
    {
        RETURN_ERROR("ERROR allocating rec_c memory", FUNC_NAME, FAILURE);
    }

    for(i_col = 0; i_col < num_samples; i_col++)
    {
        for(j = 0; j < TOTAL_IMAGE_BANDS; j++)
        {
           tmp_buf[j]  = buf[j] + i_col * num_scenes;
        }

        /*weighted fitting*/
        if (1==method)
        {
            fitting_compositing(tmp_buf, valid_datearray_scanline[i_col],
                        valid_datecount_scanline[i_col], lower_ordinal,
                        upper_ordinal, i_col, out_compositing, TRUE, TRUE, b_diagnosis, rec_c);
        }
        /*normal fitting*/
        else if (2==method)
        {
            fitting_compositing(tmp_buf, valid_datearray_scanline[i_col],
                        valid_datecount_scanline[i_col], lower_ordinal,
                        upper_ordinal, i_col, out_compositing, TRUE, FALSE, b_diagnosis, rec_c);
        }
        else if (3==method)
        {
            hot_compositing(tmp_buf, valid_datearray_scanline[i_col],
                        valid_datecount_scanline[i_col], lower_ordinal,
                        upper_ordinal, i_col, out_compositing);
        }
        else if (4==method)
        {
            average_compositing(tmp_buf, valid_datearray_scanline[i_col],
                        valid_datecount_scanline[i_col], lower_ordinal,
                        upper_ordinal, i_col, out_compositing);
        }
        else if (5==method)
        {
            fitting_compositing(tmp_buf, valid_datearray_scanline[i_col],
                        valid_datecount_scanline[i_col], lower_ordinal,
                        upper_ordinal, i_col, out_compositing, FALSE, FALSE, b_diagnosis, rec_c);
        }
        else if (6==method)
        {
            modified_hot_compositing(tmp_buf, valid_datearray_scanline[i_col],
                        valid_datecount_scanline[i_col], lower_ordinal,
                        upper_ordinal, i_col, out_compositing);
        }

        else if (7==method)
        {
            medium_compositing(tmp_buf, valid_datearray_scanline[i_col],
                               valid_datecount_scanline[i_col], lower_ordinal,
                               upper_ordinal, i_col, out_compositing);
        }
        else if (8==method)
        {
            valid_obs_count(tmp_buf, valid_datearray_scanline[i_col],
                            valid_datecount_scanline[i_col], lower_ordinal,
                            upper_ordinal, i_col, out_compositing);
        }
    }

    free(rec_c);
    free_2d_array((void **)tmp_buf);

    return SUCCESS;
}

/******************************************************************************
MODULE:  fitting_compositing

PURPOSE:  Running compositing for pixel-based time series

RETURN VALUE:
Type = int (SUCCESS, ERROR or FAILURE)

HISTORY:
Date        Programmer       Reason
--------    ---------------  -------------------------------------
05/02/2019   Su Ye         Original Development
******************************************************************************/
int fitting_compositing
(
    short int **buf,            /* I:  pixel-based time series           */
    int *valid_date_array,      /* I: valid date time series               */
    int valid_date_count,       /* I: the number of valid dates               */
    int lower_ordinal,
    int upper_ordinal,
    int i_col,
    short int **out_compositing,           /* O: outputted compositing results for four bands */
    int bfit,
    int bweighted,
    int b_diagnosis,
    Output_t* rec_c
)
{
    char FUNC_NAME[] = "fitting_compositing";
    int status;
    float date_vario;           /* I: median date                                          */
    float max_date_difference;   /* I: maximum difference between two neighbor dates        */
    float adj_rmse[TOTAL_IMAGE_BANDS]; /* Adjusted RMSE for all bands          */
    int *bl_ids;
    int n_clr;
    int n_clr_1;
    int n_outlier_1;
    int n_clr_2;
    int n_outlier_2;
    int k, b;
    int* clrx;
    float **clry;
    int* clrx_1;
    float **clry_1;
    int* clrx_2;
    float **clry_2;
    int i;
    float C0; // intercept from each test output
    float C1; // slope from each test output

    clrx = (int*)calloc(valid_date_count, sizeof(int));
    clry = (float **) allocate_2d_array (TOTAL_IMAGE_BANDS, valid_date_count,
                                             sizeof (float));
    if (clry == NULL)
    {
        RETURN_ERROR ("Allocating clry memory", FUNC_NAME, FAILURE);
    }

    /**************************************************************/
    /*                                                            */
    /*      select observations in the observation window          */
    /*                                                            */
    /**************************************************************/
    n_clr = 0;
    for(i = 0; i < valid_date_count; i++)
    {
        if((valid_date_array[i] > lower_ordinal - 1) && (valid_date_array[i] < upper_ordinal + 1))
        {
            clrx[n_clr] = valid_date_array[i];
            for(b = 0; b < TOTAL_IMAGE_BANDS; b++)
            {
                clry[b][n_clr] = (float)buf[b][i];
            }
            n_clr++;
        }
    }

    /********************************************/
    /*    condition 1: zero valid observation   */
    /********************************************/
    if(n_clr==0)
    {
        for(i = 0; i < TOTAL_IMAGE_BANDS; i++)
        {
            out_compositing[i][i_col] = -9999;
        }

        if(TRUE == b_diagnosis)
        {
            rec_c->condition = NOOBS_CONDITION;
        }

        free(clrx);
        status = free_2d_array((void **)clry);
        if (status != SUCCESS)
        {
            RETURN_ERROR ("Freeing memory: clry\n",
                          FUNC_NAME, FAILURE);
        }

        return SUCCESS;
    }

    /**********************************************/
    /*    condition 2: inefficient observations   */
    /**********************************************/
    else if (n_clr < MIN_SAMPLE)
    {
        median_compositing(buf, valid_date_array, valid_date_count,
                           i_col, out_compositing);

        if(TRUE == b_diagnosis)
        {
            rec_c->condition = INEFFICIENT_CONDITION;
        }

        free(clrx);
        status = free_2d_array((void **)clry);
        if (status != SUCCESS)
        {
            RETURN_ERROR ("Freeing memory: clry\n",
                          FUNC_NAME, FAILURE);
        }

        return SUCCESS;
    }
    else
    {
        if(TRUE == b_diagnosis)
        {
            rec_c->condition = NORMAL_CONDITION;
        }
    }

    /**********************************************/
    /*    condition 3: standard procedure         */
    /**********************************************/

    bl_ids = (int *)calloc(n_clr, sizeof(int));
    if (bl_ids == NULL)
    {
        RETURN_ERROR("ERROR allocating bl_ids memory", FUNC_NAME, FAILURE);
    }

    clrx_1 = (int *)calloc(n_clr, sizeof(int));
    clry_1 = (float **) allocate_2d_array (TOTAL_IMAGE_BANDS, n_clr,
                                             sizeof (float));
    if (clry_1 == NULL)
    {
        RETURN_ERROR ("Allocating clry_1 memory", FUNC_NAME, FAILURE);
    }

    clrx_2 = (int *)calloc(n_clr, sizeof(int));
    clry_2 = (float **) allocate_2d_array (TOTAL_IMAGE_BANDS, n_clr,
                                             sizeof (float));
    if (clry_2 == NULL)
    {
        RETURN_ERROR ("Allocating clry_2 memory", FUNC_NAME, FAILURE);
    }

    /**************************************************************/
    /*                                                            */
    /*      calculate variogram for each band and dates.          */
    /*                                                            */
    /**************************************************************/
    status = adjust_median_variogram(clrx, clry, TOTAL_IMAGE_BANDS,
                                     0, n_clr-1, &date_vario,
                                     &max_date_difference, adj_rmse);
    if (status != SUCCESS)
    {
            RETURN_ERROR("ERROR calling median_variogram routine", FUNC_NAME,
                         FAILURE);
    }

    status = greenband_test(clrx, clry, 0, n_clr-1, adj_rmse[1], T_CONST_SINGLETAIL_9999,
            bl_ids, &C0, &C1);

    if (status != SUCCESS)
    {
        RETURN_ERROR("ERROR calling greenband_test",
                      FUNC_NAME, FAILURE);
    }

    if(TRUE == b_diagnosis)
    {
        rec_c->C0_green = C0;
        rec_c->C1_green = C1;
    }

    /**************************************************/
    /*                                                */
    /*         remove outliers.                       */
    /*                                                */
    /**************************************************/
    n_clr_1 = 0;
    n_outlier_1 = 0;
    for(i = 0; i < n_clr; i++)
    {
        if(bl_ids[i] == 0)
        {
            clrx_1[n_clr_1] = clrx[i];
            for (b = 0; b < TOTAL_IMAGE_BANDS; b++)
            {
                clry_1[b][n_clr_1] = clry[b][i];
            }
            n_clr_1 = n_clr_1 + 1;
        }
        else
        {
            if(TRUE == b_diagnosis)
            {
                rec_c->outlier_dates_green[n_outlier_1] = clrx[i];
                n_outlier_1 = n_outlier_1 + 1;
            }
        }
    }

    rec_c->n_outlier_green = n_outlier_1;

    /* if n_clr_1 < MIN_SAMPLE, means that green test failed, need to reset*/
    if (n_clr_1 < MIN_SAMPLE)
    {
        n_clr_1 = 0;
        for(i = 0; i < n_clr; i++)
        {

            clrx_1[n_clr_1] = clrx[i];
            for (b = 0; b < TOTAL_IMAGE_BANDS; b++)
            {
                clry_1[b][n_clr_1] = (float)clry[b][i];
            }
            n_clr_1 = n_clr_1 + 1;

        }
        if(TRUE == b_diagnosis)
            rec_c->b_success_green = FAILURE;
    }
    else{
        if(TRUE == b_diagnosis)
            rec_c->b_success_green = SUCCESS;
    }

    /**************************************************/
    /*                                                */
    /* nir band test                                  */
    /*                                                */
    /**************************************************/

    for (k = 0; k < n_clr_1; k++)
       bl_ids[k] = 0;

    status = nirband_test(clrx_1, clry_1, 0, n_clr_1-1, adj_rmse[3], T_CONST_SINGLETAIL_9999,
            bl_ids, &C0, &C1);

    if (status != SUCCESS)
    {
        RETURN_ERROR("ERROR calling nirband_test",
                      FUNC_NAME, FAILURE);
    }

    if(TRUE == b_diagnosis)
    {
        rec_c->C0_nir = C0;
        rec_c->C1_nir = C1;
    }

    /**************************************************/
    /*                                                */
    /*         remove outliers.                       */
    /*                                                */
    /**************************************************/
    n_clr_2 = 0;
    n_outlier_2 = 0;
    for(i = 0; i < n_clr_1; i++)
    {
        if(bl_ids[i] == 0)
        {
            clrx_2[n_clr_2] = clrx_1[i];
            for (b = 0; b < TOTAL_IMAGE_BANDS; b++)
            {
                clry_2[b][n_clr_2] = clry_1[b][i];
            }
            n_clr_2 = n_clr_2 + 1;
        }
        else
        {
            if(TRUE == b_diagnosis)
            {
                rec_c->outlier_dates_nir[n_outlier_2] = clrx_1[i];
                n_outlier_2 = n_outlier_2 + 1;
            }
        }
    }

    rec_c->n_outlier_nir = n_outlier_2;

    /* if n_clr_1 < MIN_SAMPLE, means that nir test failed, reset*/
    if (n_clr_2 < MIN_SAMPLE)
    {
        for(i = 0; i < n_clr_1; i++)
        {
            n_clr_2 = 0;
            clrx_2[n_clr_2] = clrx_1[i];
            for (b = 0; b < TOTAL_IMAGE_BANDS; b++)
            {
                clry_2[b][n_clr_2] = clry_1[b][i];
            }
            n_clr_2 = n_clr_2 + 1;
        }

        if(TRUE == b_diagnosis)
            rec_c->b_success_nir = FAILURE;
    }
    else{
        if(TRUE == b_diagnosis)
            rec_c->b_success_nir = SUCCESS;
    }

    if(bfit == TRUE)
        linear_fit_centerdate(clrx_2, clry_2, n_clr_2, 0, (lower_ordinal + upper_ordinal)/2,i_col,
                          out_compositing, bweighted, rec_c->C0_final, rec_c->C1_final);
    else
    {
        float wt;
        int j;
        double index_sum[TOTAL_IMAGE_BANDS];
        for(i = 0; i < TOTAL_IMAGE_BANDS; i++)
            index_sum[i] = 0;
        float wt_sum = 0;

        for(i = 0; i < n_clr_2; i++)
        {
            wt = (float)1/((clry_2[BLUE_INDEX][i] - 0.5 * clry_2[RED_INDEX][i]) * (clry_2[BLUE_INDEX][i] - 0.5 * clry_2[RED_INDEX][i]));
            for(j = 0; j < TOTAL_IMAGE_BANDS; j++)
            {
                index_sum[j] = index_sum[j] + clry_2[j][i] * wt;
            }
            wt_sum = wt_sum + wt;
        }

        /********************************************/
        /*    condition 2: standard procedures      */
        /********************************************/

        for(j = 0; j < TOTAL_IMAGE_BANDS; j++)
        {
            out_compositing[j][i_col] = (short int)(index_sum[j] / wt_sum);
        }

    }

    free(bl_ids);
    free(clrx);
    status = free_2d_array((void **)clry);
    if (status != SUCCESS)
    {
        RETURN_ERROR ("Freeing memory: clry\n",
                      FUNC_NAME, FAILURE);
    }

    free(clrx_1);
    status = free_2d_array((void **)clry_1);
    if (status != SUCCESS)
    {
        RETURN_ERROR ("Freeing memory: clry_1\n",
                      FUNC_NAME, FAILURE);
    }

    free(clrx_2);
    status = free_2d_array((void **)clry_2);
    if (status != SUCCESS)
    {
        RETURN_ERROR ("Freeing memory: clry_2\n",
                      FUNC_NAME, FAILURE);
    }

    return SUCCESS;

}

//int fitting_compositing_scanline
//(
//    short int **buf,            /* I:  scanline-based time series           */
//    int **valid_datearray_scanline,      /* I: valid date time series               */
//    int *valid_datecount_scanline,       /* I: the number of valid dates               */
//    int lower_ordinal,                /* I: lower ordinal dates               */
//    int upper_ordinal,                   /* I: upper_ordinal for temporal range of composition   */
//    int num_samples,                /* I: the pixel number in a row              */
//    int num_scenes,                 /* I: the number of scenes               */
//    short int **out_compositing           /* O: outputted compositing results for four bands */
//)
//{
//    int i, j;
//    int i_col;
//    short int **tmp_buf;                   /* This is the image bands buffer, valid pixel only*/
//    int result;
//    char FUNC_NAME[] = "compositing_scanline";
//    tmp_buf = (short int **) allocate_2d_array (TOTAL_IMAGE_BANDS, num_scenes, sizeof (short int));
//    if(tmp_buf == NULL)
//    {
//        RETURN_ERROR("ERROR allocating tmp_buf memory", FUNC_NAME, FAILURE);
//    }

//    for(i_col = 0; i_col < num_samples; i_col++)
//    {
//        for(j = 0; j < TOTAL_IMAGE_BANDS; j++)
//        {
//           tmp_buf[j]  = buf[j] + i_col * num_scenes;
//        }

//        fitting_compositing(tmp_buf, valid_datearray_scanline[i_col],
//                    valid_datecount_scanline[i_col], lower_ordinal,
//                    upper_ordinal, i_col, out_compositing);
//    }

//    free_2d_array((void **)tmp_buf);

//    return SUCCESS;
//}

ext: c

max_stars_repo_name: agroimpacts/imager

max_issues_repo_name: agroimpacts/imager

max_forks_repo_name: agroimpacts/imager


text:
////////////////////////////////////////////////////////////
//
// Copyright (c) 2018 Jan Filipowicz, Filip Turobos
//
// 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 SALESMAN_EXAMPLE_PATH_EVALUATOR_H
#define SALESMAN_EXAMPLE_PATH_EVALUATOR_H

#include <functional>
#include <numeric>
#include <gsl/gsl_assert>
#include "permutation.h"

template<class Matrix>
class path_evaluator {
public:
	using matrix_type = Matrix;
	using row_type = typename matrix_type::value_type;
	using value_type = typename row_type::value_type;
	explicit path_evaluator(const matrix_type& matrix);
	value_type operator()(const permutation& perm) const;
private:
	matrix_type matrix;
};

template<class Matrix>
inline path_evaluator<Matrix>::path_evaluator(const matrix_type& matrix)
	: matrix(matrix) {}

template<class Matrix>
inline auto path_evaluator<Matrix>::operator()(const permutation& perm) const -> value_type {
	Expects(perm.size() == matrix.size());
	Expects(perm.size() > 0);
	const auto distance = [this](unsigned dest, unsigned src) {
		return gsl::at(gsl::at(matrix, src), dest);
	};
	return std::inner_product(std::next(perm.begin()), perm.end(), perm.begin(), distance(perm.front(), perm.back()), std::plus<>(), distance);
}

#endif

max_stars_repo_name: SirEmentaler/Eugenics-Wars

max_issues_repo_name: SirEmentaler/Eugenics-Wars

max_forks_repo_name: SirEmentaler/Eugenics-Wars


text:
/* usleep, inet_aton */
#define _GNU_SOURCE

#include <stdio.h>
#include <time.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
#include <pthread.h>
#include <fcntl.h>
#include <assert.h>
#include <sys/types.h>
#include <sys/socket.h>
#include <netinet/in.h>
#include <arpa/inet.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include "agent.h"

/*
 * Grafo con 6 nodi -> dmax = 5
 * devo prendere una frequenza  fs > 2*dmax/pi
 * e quindi una ts < pi/(2*dmax) = 0.314 sec
 * prendiamo h di runge kutta = ts
 */

int _agents_number = 6;

float h = 0.05;

int N = 32768;

double x0[6] = {0.00, 0.00, 0.00, 0.00, 0.00, 0.00};

double z0[6] = {1.00, 0.00, 0.00, 0.00, 0.00, 0.00};

typedef enum {
    log_debug = 1,
    log_normal = 2,
    log_details = 4
} log_t;

int log_level = log_normal;

/*int log_level = log_debug | log_normal | log_details;*/

struct msg_t {
    int msg_id;
    int msg_time;
    float msg_x;
    float msg_z;
};

struct agent_t {
    int agent_id;
    int agent_x;
    int agent_z;
    pthread_mutex_t agent_mutex;
};

typedef struct msg_t *msg;

static pthread_mutex_t agent_class_mutex;
static int agentObjectsNumber = 0;

static void initClassAgent() {
    pthread_mutex_init(&(agent_class_mutex), NULL);
}

static void cleanClassAgent() {
    pthread_mutex_destroy(&(agent_class_mutex));

}

agent allocAgent() {

    agent agent_M_;

    if (agentObjectsNumber == 0) initClassAgent();

    agentObjectsNumber++;

    agent_M_ = (agent) malloc(sizeof (struct agent_t));

    agent_M_->agent_id = 0;
    pthread_mutex_init(&(agent_M_->agent_mutex), NULL);

    return agent_M_;

}

void freeAgent(agent _F_agent) {

    pthread_mutex_destroy(&(_F_agent->agent_mutex));

    free(_F_agent);

    _F_agent = NULL;

    agentObjectsNumber--;

    if (agentObjectsNumber == 0) cleanClassAgent();
}

void setIdAgent(agent agent_p) {
    agent_p = NULL;
}

void getIdAgent(agent agent_p) {
    agent_p = NULL;
}

void setAgentState(agent agent_p) {
    agent_p = NULL;
}

void getAgentState(agent agent_p) {
    agent_p = NULL;
}

static int connectTo(char *_ip_address_remote, int _ip_port_remote) {

    int socket_M_;

    struct sockaddr_in sockaddr_in_remote_;

    memset(&sockaddr_in_remote_, 0, sizeof (sockaddr_in_remote_));
    sockaddr_in_remote_.sin_family = AF_INET;
    sockaddr_in_remote_.sin_port = htons(_ip_port_remote);
    inet_aton(_ip_address_remote, &(sockaddr_in_remote_.sin_addr));

    if ((socket_M_ = socket(AF_INET, SOCK_STREAM, 0)) == -1) {
        perror("connectTo: socket() Error\n");
        exit(1);
    }

    while (connect(socket_M_, (struct sockaddr *) & sockaddr_in_remote_, sizeof (struct sockaddr_in)) == -1) {
        perror("connectTo: connect() Error\n");
        exit(1);
    }

    return socket_M_;
}

static void sendMsg(int _socket, msg _msg) {

    if (send(_socket, _msg, sizeof (struct msg_t), 0) == -1) {
        perror("sendMsg: send() Error\n");
        exit(1);
    };
}

static msg recvMsg(int _socket) {

    msg _M_msg = (struct msg_t *) malloc(sizeof (struct msg_t));

    if (recv(_socket, _M_msg, sizeof (struct msg_t), 0) == -1) {
        perror("recvMsg: recv() Error\n");
        exit(1);
    }

    return _M_msg;
}

parametri_agent allocParamAgent(int _identity, char *_agent_router_ip, int _port_agent_to_router, int _port_router_to_agent) {

    parametri_agent _M_parametri_agent = (parametri_agent) malloc(sizeof (struct parametri_agent_t));

    _M_parametri_agent->identity = _identity;
    _M_parametri_agent->agent_router_ip = _agent_router_ip;
    _M_parametri_agent->port_router_to_agent = _port_router_to_agent;
    _M_parametri_agent->port_agent_to_router = _port_agent_to_router;

    return _M_parametri_agent;
}

void *runAgent(void *_parametri_agent) {

    parametri_agent _parametri = (parametri_agent) _parametri_agent;

    int _identity = _parametri->identity;

    char *_agent_router_ip;
    int _port_router_to_agent;
    int _port_agent_to_router;

    int agent_to_router_socket_;
    int router_to_agent_socket_;

    int _agent_neighborhood;

    int i, c;

    msg _M_msg_in, _M_msg_out;

    double x[4] = {0.00, 0.00, 0.00, 0.00};
    double z[4] = {0.00, 0.00, 0.00, 0.00};
    double r[4] = {0.00, 0.00, 0.00, 0.00};
    double sx[4] = {0.00, 0.00, 0.00, 0.00};
    double sz[4] = {0.00, 0.00, 0.00, 0.00};

    FILE *agent_log;
    char agent_log_file_name[FILENAME_MAX];

    gsl_matrix *_adjacency_matrix;
    FILE *agent_adjacency_matrix_file;

    _adjacency_matrix = gsl_matrix_alloc(_agents_number, _agents_number);
    agent_adjacency_matrix_file = fopen("adjacency_matrix.txt", "r");
    gsl_matrix_fscanf(agent_adjacency_matrix_file, _adjacency_matrix);
    fclose(agent_adjacency_matrix_file);

    /*
     * Apertura del file di log
     */

    sprintf(agent_log_file_name, "agent_%d.log", _identity);

    agent_log = fopen(agent_log_file_name, "w");

    if (agent_log == NULL) {
        perror("Errore apertura file: \"agent.log\"\n");
        exit(1);
    };

    setbuf(agent_log, NULL);

    if (log_level & log_debug) fprintf(agent_log, "A%d: Open log\n", _identity);

    if (log_level & log_debug) fprintf(agent_log, "A%d: Parametri di invocazione:\n", _identity);

    if (log_level & log_debug) fprintf(agent_log, "A%d: > _identity = %d\n", _identity, _identity);

    _agent_router_ip = _parametri->agent_router_ip;

    if (log_level & log_debug) fprintf(agent_log, "A%d: > _agent_router_ip = %s\n", _identity, _agent_router_ip);

    _port_agent_to_router = _parametri->port_agent_to_router;

    if (log_level & log_debug) fprintf(agent_log, "A%d: > _port_agent_to_router = %d\n", _identity, _port_agent_to_router);

    _port_router_to_agent = _parametri->port_router_to_agent;

    if (log_level & log_debug) fprintf(agent_log, "A%d: > _port_router_to_agent = %d\n", _identity, _port_router_to_agent);

    if (log_level & log_debug) fprintf(agent_log, "A%d: Memoria:\n", _identity);

    if (log_level & log_debug) fprintf(agent_log, "A%d: > h = %2.8g\n", _identity, h);

    /*
     * Calcolo del numero di agent vicini attraverso l'analisi della adjacency_matrix
     */

    _agent_neighborhood = 0;

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

        _agent_neighborhood = _agent_neighborhood + (int) gsl_matrix_get(_adjacency_matrix, i, _identity - 1);

        if (log_level & log_debug) fprintf(agent_log, "A%d: > gsl_matrix_get(%d,%d) = %d\n", _identity, i, _identity - 1, (int) gsl_matrix_get(_adjacency_matrix, i, _identity - 1));
    };

    if (log_level & log_debug) fprintf(agent_log, "A%d: > _agent_neighborhood = %d\n", _identity, _agent_neighborhood);

    if (log_level & log_debug) fprintf(agent_log, "A%d: > _agents_number = %d\n", _identity, _agents_number);

    _M_msg_out = (struct msg_t *) malloc(sizeof (struct msg_t));

    if (log_level & log_debug) fprintf(agent_log, "A%d: > alloc(msg_out)\n", _identity);

    if (log_level & log_debug) fprintf(agent_log, "A%d: Socket:\n", _identity);

    agent_to_router_socket_ = connectTo(_agent_router_ip, _port_agent_to_router);

    if (log_level & log_debug) fprintf(agent_log, "A%d: > connectTo(%s,%d)\n", _identity, _agent_router_ip, _port_agent_to_router);

    router_to_agent_socket_ = connectTo(_agent_router_ip, _port_router_to_agent);

    if (log_level & log_debug) fprintf(agent_log, "A%d: > connectTo(%s,%d)\n", _identity, _agent_router_ip, _port_router_to_agent);

    if (log_level & log_debug) fprintf(agent_log, "A%d: Signal:\n", _identity);

    /*
     * Imposto l'identità dei messaggi in output
     */

    _M_msg_out->msg_id = _identity;

    /*
     * Imposto il timer interno
     */

    _M_msg_out->msg_time = 0;

    /*
     * Invio di un messaggio contenente la propria identità sui due canali di in/out
     */

    _M_msg_out->msg_x = 0.00;
    _M_msg_out->msg_z = 0.00;

    sendMsg(router_to_agent_socket_, _M_msg_out);

    if (log_level & log_debug) fprintf(agent_log, "A%d: > snd (%d,%d,%2.8g,%2.8g) > %d\n", _identity, _M_msg_out->msg_id, _M_msg_out->msg_time, _M_msg_out->msg_x, _M_msg_out->msg_z, _port_router_to_agent);

    sendMsg(agent_to_router_socket_, _M_msg_out);

    if (log_level & log_debug) fprintf(agent_log, "A%d: > snd (%d,%d,%2.8g,%2.8g) > %d\n", _identity, _M_msg_out->msg_id, _M_msg_out->msg_time, _M_msg_out->msg_x, _M_msg_out->msg_z, _port_agent_to_router);

    /*
     * Ciclo principale
     */

    if (log_level & log_debug) fprintf(agent_log, "A%d: Loop:\n", _identity);

    if (log_level & log_debug) fprintf(agent_log, "A%d: > for i=[0,%d]\n", _identity, N);

    if (log_level & log_debug) fprintf(agent_log, "A%d: ----- Loop Start -----\n", _identity);

    x[0] = x0[_identity - 1];
    z[0] = z0[_identity - 1];

    for (c = 0; c < N; c++) {

        if (log_level & log_debug) fprintf(agent_log, "A%d: > i = %d\n", _identity, c);

        /*
         * snd x
         */

        _M_msg_out->msg_time = 0;
        _M_msg_out->msg_x = x[0];
        _M_msg_out->msg_z = z[0];
        sendMsg(agent_to_router_socket_, _M_msg_out);
        if (log_level & log_normal) fprintf(agent_log, "%2.8g %2.8g\n", _M_msg_out->msg_x, _M_msg_out->msg_z);
        if (log_level & log_details) fprintf(agent_log, "A%d: > snd (%d,%d,%2.8g,%2.8g)\n", _identity, _M_msg_out->msg_id, _M_msg_out->msg_time, _M_msg_out->msg_x, _M_msg_out->msg_z);

        /*
         * k1
         */

        while (r[0] < _agent_neighborhood) {
            _M_msg_in = recvMsg(router_to_agent_socket_);
            if (log_level & log_details) fprintf(agent_log, "A%d: > rcv (%d,%d,%2.8g,%2.8g)\n", _identity, _M_msg_in->msg_id, _M_msg_out->msg_time, _M_msg_in->msg_x, _M_msg_in->msg_z);
            sx[_M_msg_in->msg_time] = sx[_M_msg_in->msg_time] + _M_msg_in->msg_x;
            sz[_M_msg_in->msg_time] = sz[_M_msg_in->msg_time] + _M_msg_in->msg_z;
            r[_M_msg_in->msg_time]++;
        }

        x[1] = r[0] * z[0] - sz[0];
        z[1] = sx[0] - r[0] * x[0];
        sx[0] = 0;
        sz[0] = 0;
        r[0] = 0;

        _M_msg_out->msg_time = 1;
        _M_msg_out->msg_x = x[1];
        _M_msg_out->msg_z = z[1];
        sendMsg(agent_to_router_socket_, _M_msg_out);
        if (log_level & log_details) fprintf(agent_log, "A%d: > snd (%d,%d,%2.8g,%2.8g)\n", _identity, _M_msg_out->msg_id, _M_msg_out->msg_time, _M_msg_out->msg_x, _M_msg_out->msg_z);

        /*
         * k2
         */

        while (r[1] < _agent_neighborhood) {
            _M_msg_in = recvMsg(router_to_agent_socket_);
            if (log_level & log_details) fprintf(agent_log, "A%d: > rcv (%d,%d,%2.8g,%2.8g)\n", _identity, _M_msg_in->msg_id, _M_msg_out->msg_time, _M_msg_in->msg_x, _M_msg_in->msg_z);
            sx[_M_msg_in->msg_time] = sx[_M_msg_in->msg_time] + _M_msg_in->msg_x;
            sz[_M_msg_in->msg_time] = sz[_M_msg_in->msg_time] + _M_msg_in->msg_z;
            r[_M_msg_in->msg_time]++;
        }

        x[2] = x[1] + h * (r[1] * z[1] - sz[1]) / 2;
        z[2] = z[1] + h * (sx[1] - r[1] * x[1]) / 2;
        sx[1] = 0;
        sz[1] = 0;
        r[1] = 0;

        _M_msg_out->msg_time = 2;
        _M_msg_out->msg_x = x[2];
        _M_msg_out->msg_z = z[2];
        sendMsg(agent_to_router_socket_, _M_msg_out);
        if (log_level & log_details) fprintf(agent_log, "A%d: > snd (%d,%d,%2.8g,%2.8g)\n", _identity, _M_msg_out->msg_id, _M_msg_out->msg_time, _M_msg_out->msg_x, _M_msg_out->msg_z);

        /*
         * k3
         */

        while (r[2] < _agent_neighborhood) {
            _M_msg_in = recvMsg(router_to_agent_socket_);
            if (log_level & log_details) fprintf(agent_log, "A%d: > rcv (%d,%d,%2.8g,%2.8g)\n", _identity, _M_msg_in->msg_id, _M_msg_out->msg_time, _M_msg_in->msg_x, _M_msg_in->msg_z);
            sx[_M_msg_in->msg_time] = sx[_M_msg_in->msg_time] + _M_msg_in->msg_x;
            sz[_M_msg_in->msg_time] = sz[_M_msg_in->msg_time] + _M_msg_in->msg_z;
            r[_M_msg_in->msg_time]++;
        }

        x[3] = x[1] + h * (r[2] * z[2] - sz[2]) / 2;
        z[3] = z[1] + h * (sx[2] - r[2] * x[2]) / 2;
        sx[2] = 0;
        sz[2] = 0;
        r[2] = 0;

        _M_msg_out->msg_time = 3;
        _M_msg_out->msg_x = x[3];
        _M_msg_out->msg_z = z[3];
        sendMsg(agent_to_router_socket_, _M_msg_out);
        if (log_level & log_details) fprintf(agent_log, "A%d: > snd (%d,%d,%2.8g,%2.8g)\n", _identity, _M_msg_out->msg_id, _M_msg_out->msg_time, _M_msg_out->msg_x, _M_msg_out->msg_z);

        /*
         * k4, x
         */

        while (r[3] < _agent_neighborhood) {
            _M_msg_in = recvMsg(router_to_agent_socket_);
            if (log_level & log_details) fprintf(agent_log, "A%d: > rcv (%d,%d,%2.8g,%2.8g)\n", _identity, _M_msg_in->msg_id, _M_msg_out->msg_time, _M_msg_in->msg_x, _M_msg_in->msg_z);
            sx[_M_msg_in->msg_time] = sx[_M_msg_in->msg_time] + _M_msg_in->msg_x;
            sz[_M_msg_in->msg_time] = sz[_M_msg_in->msg_time] + _M_msg_in->msg_z;
            r[_M_msg_in->msg_time]++;
        }

        x[0] = x[0] + h * (x[1] + x[2] + x[3] + h * (r[3] * z[3] - sz[3]) / 2) / 3;
        z[0] = z[0] + h * (z[1] + z[2] + z[3] + h * (sx[3] - r[3] * x[3]) / 2) / 3;
        sx[3] = 0;
        sz[3] = 0;
        r[3] = 0;

    }

    if (log_level & log_debug) fprintf(agent_log, "A%d: ----- Loop End ------\n", _identity);

    if (log_level & log_debug) fprintf(agent_log, "A%d: Signal:\n", _identity);

    /*
     * Invio del killer message
     */

    _M_msg_out->msg_id = -1;
    _M_msg_out->msg_x = 0;
    _M_msg_out->msg_z = 0;

    sendMsg(agent_to_router_socket_, _M_msg_out);

    if (log_level & log_debug) fprintf(agent_log, "A%d: > snd (%d,%d,%2.8g,%2.8g) > %d\n", _identity, _M_msg_out->msg_id, _M_msg_out->msg_time, _M_msg_out->msg_x, _M_msg_out->msg_z, _port_agent_to_router);

    if (log_level & log_debug) fprintf(agent_log, "A%d: Socket:\n", _identity);

    /*
     * Chiusura delle sockets
     */

    close(router_to_agent_socket_);

    if (log_level & log_debug) fprintf(agent_log, "A%d: > close(%d)\n", _identity, _port_router_to_agent);

    close(agent_to_router_socket_);

    if (log_level & log_debug) fprintf(agent_log, "A%d: > close(%d)\n", _identity, _port_agent_to_router);

    if (log_level & log_debug) fprintf(agent_log, "A%d: Memoria:\n", _identity);

    /*
     * Free delle variabili usate
     */

    free(_M_msg_out);

    if (log_level & log_debug) fprintf(agent_log, "A%d: > free(msg_out)\n", _identity);

    free(_parametri_agent);

    if (log_level & log_debug) fprintf(agent_log, "A%d: > parametri\n", _identity);

    /*
     * fflush e chiusura del logfile
     */

    if (log_level & log_debug) fprintf(agent_log, "A%d: Close log\n", _identity);

    fflush(agent_log);

    fclose(agent_log);

    return NULL;
}

max_stars_repo_name: vittorioc/VirtualAgent

max_issues_repo_name: vittorioc/VirtualAgent

max_forks_repo_name: vittorioc/VirtualAgent


text:
/*
Ballistic: a software to benchmam ballistic models.

AUTHORS: Javier Burguete Tolosa.

Copyright 2018, AUTHORS.

Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:

  1. Redistributions of source code must retain the above copyright notice,
    this list of conditions and the following disclaimer.

  2. Redistributions in binary form must reproduce the above copyright notice,
    this list of conditions and the following disclaimer in the
    documentation and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY AUTHORS ``AS IS'' AND ANY EXPRESS OR IMPLIED
WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
SHALL AUTHORS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
OF SUCH DAMAGE.
*/

/**
 * \file runge-kutta.c
 * \brief Source file to define the numerical method data and functions.
 * \author Javier Burguete Tolosa.
 * \copyright Copyright 2018.
 */
#define _GNU_SOURCE
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <libxml/parser.h>
#include <glib.h>
#include "config.h"
#include "utils.h"
#include "equation.h"
#include "method.h"

#define DEBUG_METHOD 0
///< macro to debug the numerical method functions.

/**
 * Function to init the numerical method.
 */
void
method_init (Method * m,        ///< Method struct.
             unsigned int nsteps,       ///< number of steps.
             unsigned int order)        ///< order of accuracy.
{
#if DEBUG_METHOD
  fprintf (stderr, "method_init: start\n");
#endif
  m->nsteps = nsteps;
  m->order = order;
#if DEBUG_METHOD
  fprintf (stderr, "method_init: end\n");
#endif
}

/**
 * Function to init the variables of the numerical method.
 */
void
method_init_variables (Method * m)      ///< Method struct.
{
  unsigned int i, n;
#if DEBUG_METHOD
  fprintf (stderr, "method_init_variables: start\n");
  fprintf (stderr, "method_init_variables: nsteps=%u\n", m->nsteps);
#endif
  n = m->nsteps + 1;
  m->r0 = (long double **) g_slice_alloc (n * sizeof (long double *));
  m->r1 = (long double **) g_slice_alloc (n * sizeof (long double *));
  m->r2 = (long double **) g_slice_alloc (n * sizeof (long double *));
  for (i = 0; i < n; ++i)
    {
      m->r0[i] = (long double *) g_slice_alloc (3 * sizeof (long double));
      m->r1[i] = (long double *) g_slice_alloc (3 * sizeof (long double));
      m->r2[i] = (long double *) g_slice_alloc (3 * sizeof (long double));
    }
  m->et0 = m->et1 = 0.L;
#if DEBUG_METHOD
  fprintf (stderr, "method_init_variables: end\n");
#endif
}

/**
 * Function to calculate the following numerical step size based on error
 * control.
 *
 * \return next time step size.
 */
long double
method_dt (Method * m,          ///< Method struct.
           long double dt)      ///< actual time step size.
{
  long double dt2;
#if DEBUG_METHOD
  fprintf (stderr, "method_dt: start\n");
#endif
  dt2 = dt * fminl (m->alpha,
                    powl (m->emt * dt / m->e0, 1.L / (m->order - 1.L)));
#if DEBUG_METHOD
  fprintf (stderr, "method_dt: dt=%Lg\n", dt2);
  fprintf (stderr, "method_dt: end\n");
#endif
  return dt2;
}

/**
 * Function to free the memory used by the Method struct.
 */
void
method_delete (Method * m)      ///< Method struct.
{
  unsigned int i, n;
#if DEBUG_METHOD
  fprintf (stderr, "method_delete: start\n");
  fprintf (stderr, "method_delete: nsteps=%u\n", m->nsteps);
#endif
  n = i = m->nsteps + 1;
  do
    {
      --i;
      g_slice_free1 (3 * sizeof (long double), m->r2[i]);
      g_slice_free1 (3 * sizeof (long double), m->r1[i]);
      g_slice_free1 (3 * sizeof (long double), m->r0[i]);
    }
  while (i);
  g_slice_free1 (n * sizeof (long double *), m->r2);
  g_slice_free1 (n * sizeof (long double *), m->r1);
  g_slice_free1 (n * sizeof (long double *), m->r0);
#if DEBUG_METHOD
  fprintf (stderr, "method_delete: end\n");
#endif
}

/**
 * Function to read the numerical method data on a XML node.
 *
 * \return 1 on success, 0 on error.
 */
int
method_read_xml (Method * m,    ///< Method struct.
                 xmlNode * node)        ///< XML node.
{
  const char *message[] = {
    "Bad dt",
    "Bad alpha",
    "Bad beta",
    "Bad error per time",
    "Unknown error control type"
  };
  int e, error_code;
#if DEBUG_METHOD
  fprintf (stderr, "method_read_xml: start\n");
#endif
  m->error_dt = xml_node_get_uint (node, XML_TIME_STEP, &error_code);
  if (error_code)
    {
      e = 0;
      goto fail;
    }
  switch (m->error_dt)
    {
    case 0:
      m->emt = 0.L;
      break;
    case 1:
      m->alpha = xml_node_get_float (node, XML_ALPHA, &error_code);
      if (error_code)
        {
          e = 1;
          goto fail;
        }
      m->beta = xml_node_get_float (node, XML_BETA, &error_code);
      if (error_code)
        {
          e = 2;
          goto fail;
        }
      m->emt = xml_node_get_float (node, XML_ERROR_TIME, &error_code);
      if (error_code)
        {
          e = 3;
          goto fail;
        }
      break;
    default:
      e = 4;
      goto fail;
    }
#if DEBUG_METHOD
  fprintf (stderr, "method_read_xml: success\n");
  fprintf (stderr, "method_read_xml: end\n");
#endif
  return 1;

fail:
  error_add (message[e]);
#if DEBUG_METHOD
  fprintf (stderr, "method_read_xml: error\n");
  fprintf (stderr, "method_read_xml: end\n");
#endif
  return 0;
}

max_stars_repo_name: jburguete/ballistic

max_issues_repo_name: jburguete/ballistic

max_forks_repo_name: jburguete/ballistic


text:
/*
 * File: spce_sect.h
 * External definitions for spce_sect.c
 *
 *
 */

#ifndef spce_sect_H
#define spce_sect_H

#include <math.h>
#include <gsl/gsl_roots.h>
#include "aXe_grism.h"
#include "spc_trace_functions.h"

/**
  @package spce_sect
*/

/**
    A structure to support old gsl 0.9 gsl_interval
*/
typedef struct
{
     double lower;
     double upper;
}
gsl_interval;

/**
  A structure describing the function that has a zero at the section point,
  including the measuring point and the slope of the line through it,
  the function describing the spectrum trace, functions related to
  gsl and a flag to initiate special handling of vertically oriented
  objects.

  @see fill_in_sectionfun
  @see free_sectionfun
*/
typedef struct
{
  int vertical;	     /* Special handling if orientation is close to vertical */
  double m;	     /* slope of line through x0, y0 that is to */
  double x0, y0;     /* intersect the trace */
  trace_func *func;  /* Parametrization of the trace */

     /* the GSL stuff has to be kept in here to avoid excessive re-allocing
        of the solver for each pixel */
  gsl_interval *interv;
  gsl_function *gslfun;
  gsl_root_fsolver *solver;
}
sectionfun;

/* public */

int fill_in_sectionfun (sectionfun * const sf, const double inclination,
			const beam * const b);
int find_section_point (sectionfun *sf, const double x, const double y,
			double *const res);
void free_sectionfun (sectionfun * const sf);
#endif

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


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

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_exp.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_bessel.h>

#include "error.h"

#include "bessel.h"
#include "bessel_temme.h"

/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/

int
gsl_sf_bessel_Inu_scaled_e(double nu, double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x < 0.0 || nu < 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x*x < 10.0*(nu+1.0)) {
    gsl_sf_result b;
    double ex = exp(-x);
    int stat = gsl_sf_bessel_IJ_taylor_e(nu, x, 1, 100, GSL_DBL_EPSILON, &b);
    result->val  = b.val * ex;
    result->err  = b.err * ex;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return stat;
  }
  else if(0.5/(nu*nu + x*x) < GSL_ROOT3_DBL_EPSILON) {
    return gsl_sf_bessel_Inu_scaled_asymp_unif_e(nu, x, result);
  }
  else {
    int N = (int)(nu + 0.5);
    double mu = nu - N;      /* -1/2 <= mu <= 1/2 */
    double K_mu, K_mup1, Kp_mu;
    double K_nu, K_nup1, K_num1;
    double I_nu_ratio;
    int stat_Irat;
    int stat_Kmu;
    int n;

    /* obtain K_mu, K_mup1 */
    if(x < 2.0) {
      stat_Kmu = gsl_sf_bessel_K_scaled_temme(mu, x, &K_mu, &K_mup1, &Kp_mu);
    }
    else {
      stat_Kmu = gsl_sf_bessel_K_scaled_steed_temme_CF2(mu, x, &K_mu, &K_mup1, &Kp_mu);
    }

    /* recurse forward to obtain K_num1, K_nu */
    K_nu   = K_mu;
    K_nup1 = K_mup1;

    for(n=0; n<N; n++) {
      K_num1 = K_nu;
      K_nu   = K_nup1;
      K_nup1 = 2.0*(mu+n+1)/x * K_nu + K_num1;
    }

    /* calculate I_{nu+1}/I_nu */
    stat_Irat = gsl_sf_bessel_I_CF1_ser(nu, x, &I_nu_ratio);

    /* solve for I_nu */
    result->val = 1.0/(x * (K_nup1 + I_nu_ratio * K_nu));
    result->err = GSL_DBL_EPSILON * (0.5*N + 2.0) * fabs(result->val);

    return GSL_ERROR_SELECT_2(stat_Kmu, stat_Irat);
  }
}

int
gsl_sf_bessel_Inu_e(double nu, double x, gsl_sf_result * result)
{
  gsl_sf_result b;
  int stat_I = gsl_sf_bessel_Inu_scaled_e(nu, x, &b);
  int stat_e = gsl_sf_exp_mult_err_e(x, fabs(x*GSL_DBL_EPSILON),
                                        b.val, b.err,
                                        result);
  return GSL_ERROR_SELECT_2(stat_e, stat_I);
}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_bessel_Inu_scaled(double nu, double x)
{
  EVAL_RESULT(gsl_sf_bessel_Inu_scaled_e(nu, x, &result));
}

double gsl_sf_bessel_Inu(double nu, double x)
{
  EVAL_RESULT(gsl_sf_bessel_Inu_e(nu, x, &result));
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
/** Function for dealiasing weather radar radial velocities.
 * @file libdealias.c
 * @author Adriaan Dokter, adapted from dealias.c RAVE functions by Gunther Haase
 * @date 2016-09-22
 */

#include "libdealias.h"
#include <stdio.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_vector.h>

void printDealias(const float *points, const int nDims, const float nyquist[],
    const float vradObs[], float vradDealias[], const int nPoints, const int iProfileType, const int iLayer, const int iPass){
    fprintf(stderr,"#iProfile iLayer iPass azim elev nyquist vrad vradd\n");
    for(int i=0; i<nPoints;i++){
        fprintf(stderr,"%i %i %i %3.1f %3.1f %3.1f %3.1f %3.1f\n",
        iProfileType,iLayer,iPass,points[i*nDims],points[i*nDims+1],nyquist[i],vradObs[i],vradDealias[i]);
    }
}

double test_field(float u, float v,const float *points, const float *pointsTrigon, const int nPoints, const int nDims, double x[], double y[], const float nyquist[]){
    double xt, yt, e, vm;
    double esum = 0;
    for (int iPoint=0; iPoint<nPoints; iPoint++) {
        // calculate the radial velocities for the test wind fields, eq 4 in Haase et al. 2004 jaot
        vm = (u*pointsTrigon[3*iPoint] + v*pointsTrigon[3*iPoint+1])*pointsTrigon[3*iPoint+2];
        // Eq. 7 for the test radial wind:
        xt = nyquist[iPoint]/M_PI * cos(vm*M_PI/nyquist[iPoint]);
        // Eq. 6 for the test radial wind:
        yt = nyquist[iPoint]/M_PI * sin(vm*M_PI/nyquist[iPoint]);

        // summed absolute differences between observed velocity field and test fields
        e = fabs(xt-x[iPoint]) + fabs(yt-y[iPoint]);
//        fprintf(stderr,"test_field: i=%i,u=%f,v=%f,vm=%f,xt=%f,yt=%f,e=%f,sin=%f,cos=%f,cos=%f\n",iPoint,u,v,vm,xt,yt,e,pointsTrigon[3*iPoint],pointsTrigon[3*iPoint+1],pointsTrigon[3*iPoint+2]);
        if (!isnan(e)) {
            esum = esum + e;
        }
    }
    return esum;
}

double test_field_gsl(const gsl_vector *uv, void* params){
    double u,v;
    float *points = ((void **) params)[0];
    float *pointsTrigon = ((void **) params)[1];
    int nPoints = *(int *) ((void **) params)[2];
    int nDims = *(int *) ((void **) params)[3];
    double *x = ((void **) params)[4];
    double *y = ((void **) params)[5];
    float *nyquist = ((void **) params)[6];

    u = gsl_vector_get(uv, 0);
    v = gsl_vector_get(uv, 1);

    return test_field(u, v, points, pointsTrigon, nPoints, nDims, x, y, nyquist);

}

int fit_field_gsl(gsl_vector *uv, void *params){
    gsl_vector *ss;

    // define which minimizer to use
    const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex2;
    gsl_multimin_fminimizer *s = NULL;

    int iter = 0;
    int status;
    double size;
    double u1 = 0;
    double v1 = 0;

    // Set initial step sizes to 1
    ss = gsl_vector_alloc(2);
    gsl_vector_set_all(ss, 1);

    // Initialize method
    gsl_multimin_function minex_func;
    minex_func.n = 2;
    minex_func.f = &test_field_gsl;
    minex_func.params = params;
    s = gsl_multimin_fminimizer_alloc (T, 2);
    gsl_multimin_fminimizer_set (s, &minex_func, uv, ss);

    // minimize by iteration
    do
    {
        iter++;
        status = gsl_multimin_fminimizer_iterate(s);

        if (status)
        break;

        size = gsl_multimin_fminimizer_size (s);
        status = gsl_multimin_test_size (size, 1e-2);

        if (status == GSL_SUCCESS)
        {
            #ifdef FPRINTFON
            printf ("converged to minimum on iteration ");
            printf ("%d at %10.3e %10.3e f() = %7.3f size = %.3f\n",
                iter,
                gsl_vector_get (s->x, 0),
                gsl_vector_get (s->x, 1),
                s->fval, size);
            #endif

            u1=gsl_vector_get(s->x, 0);
            v1=gsl_vector_get(s->x, 1);

            gsl_vector_set(uv, 0, u1);
            gsl_vector_set(uv, 1, v1);
        }
    }
    while (status == GSL_CONTINUE && iter < 100);

    #ifdef FPRINTFON
    fprintf(stdout,"Finished dealias at (x,y)=%f,%f at f()=%f ...\n",u1,v1,s->fval);
    #endif

    // clean up
    gsl_vector_free(ss);
    gsl_multimin_fminimizer_free (s);

    if (status != GSL_SUCCESS) return 0;
    return 1;
}

int dealias_points(const float *points, const int nDims, const float nyquist[],
    const double NI_MIN, const float vo[], float vradDealias[], const int nPoints){

    int i, j, n, m, eind, fitOk;
    double min1, esum, u1, v1, min2, dmy;

    // number of rows
    m = DEALIAS_VAF;
    // number of columns
    n = DEALIAS_NF;

    // max number of folds of nyquist interval to test for
    double MVA=2*ceil(DEALIAS_VMAX/(2*NI_MIN));

    // polarscan matrix, torus projected x coordinate, eq. 6 Haase et al. 2004 jaot
    double *x = RAVE_CALLOC ((size_t)nPoints, sizeof(double));
    // polarscan matrix, torus projected y coordinate, eq. 7 Haase et al. 2004 jaot
    double *y = RAVE_CALLOC ((size_t)nPoints, sizeof(double));
    // U-components of test velocity fields
    double *uh = RAVE_CALLOC ((size_t)(m*n), sizeof(double));
    // V-components of test velocity fields
    double *vh = RAVE_CALLOC ((size_t)(m*n), sizeof(double));
    // radial velocities of the best fitting test field
    double *vt1 = RAVE_CALLOC ((size_t)nPoints, sizeof(double));
    // array with trigonometric conversions of the points array
    float *pointsTrigon = RAVE_CALLOC ((size_t)(3*nPoints), sizeof(double));

    // map measured data to 3D
    for (i=0; i<nPoints; i++) {
        x[i] = nyquist[i]/M_PI * cos(vo[i]*M_PI/nyquist[i]);
        y[i] = nyquist[i]/M_PI * sin(vo[i]*M_PI/nyquist[i]);
    }

    // trigonometric conversion points array (compute once to speed up code)
    for (int iPoint=0; iPoint<nPoints; iPoint++) {
        pointsTrigon[3*iPoint+0] = sin(points[nDims*iPoint]*DEG2RAD);
        pointsTrigon[3*iPoint+1] = cos(points[nDims*iPoint]*DEG2RAD);
        pointsTrigon[3*iPoint+2] = cos(points[nDims*iPoint+1]*DEG2RAD);
    }

    // Setting up the u and v component of the test velocity fields:
    // index n=DEALIAS_NF gives number of azimuthal directions (default n=40, i.e. steps of 360/40=9 degrees)
    // index m=DEALIAS_VAF/NI_MIN*DEALIAS_VMAX gives number of speeds (maximum speed is DEALIAS_VMAX, steps of NI_MIN/DEALIAS_VAF)
    for (i=0; i<n; i++) {
        for (j=0; j<m; j++) {
            *(uh+i*m+j) = DEALIAS_VMAX/DEALIAS_VAF*(j+1) * sin(2*M_PI/DEALIAS_NF*i);
            *(vh+i*m+j) = DEALIAS_VMAX/DEALIAS_VAF*(j+1) * cos(2*M_PI/DEALIAS_NF*i);
        }
    }

    min1 = 1e32;
    eind = 0;
    u1 = 0;
    v1 = 0;

    gsl_vector *uv;
    uv = gsl_vector_alloc(2);

    void *params[7] = {(void *) points, (void *) pointsTrigon, (void *) &nPoints, (void *) &nDims, (void *) x, (void *) y, (void *) nyquist};

    // try several test velocity fields for use as starting point in GSL fit

    for (i=0; i<m*n; i++) {

        gsl_vector_set(uv, 0, *(uh+i));
        gsl_vector_set(uv, 1, *(vh+i));
        esum = test_field_gsl(uv, &params);

        if (esum<min1) {
            min1 = esum;
            eind = i;
        }
        u1 = *(uh+eind);
        v1 = *(vh+eind);
    }
    gsl_vector_set(uv, 0, u1);
    gsl_vector_set(uv, 1, v1);

    #ifdef FPRINTFON
    fprintf(stdout,"Start dealiasing at (x,y)=%f,%f at f()=%f ...\n",u1,v1,esum);
    #endif

    fitOk = fit_field_gsl(uv, &params);
    if(!fitOk) goto cleanup;

    // the radial velocity of the best fitting test velocity field:
    for (int iPoint=0; iPoint<nPoints; iPoint++) {
        *(vt1+iPoint) = (u1*sin(points[nDims*iPoint]*DEG2RAD) + v1*cos(points[nDims*iPoint]*DEG2RAD))
                        *cos(points[nDims*iPoint+1]*DEG2RAD);
    }

    // dealias the observed velocities using the best test velocity field
    for (int iPoint=0; iPoint<nPoints; iPoint++) {
        min2 = 1e32;
        dmy = 0;
        float diffVTest = (*(vt1+iPoint)-*(vo+iPoint));
        float dv = 0;
        for (i=0; i<MVA+1; i++) {
            // get a candidate velocity fold, i.e. integer multiple of nyquist velocity
            dv = nyquist[iPoint]*(2*i-MVA);
            // checking how many folds we have, vt1-vo is the residual between the real velocity
            // and the folded velocity; which equals a  multiple of the nyquist interval
            dmy = fabs(dv-diffVTest);
            if ((dmy<min2) && (!isnan(dmy))) {
                // add the aliased interval to the observed velocity field, and obtain dealiased velocity
                *(vradDealias+iPoint) = *(vo+iPoint) + dv;
                min2 = dmy;
            }
        } // loop MVA
    } // loop over points

    cleanup:
        RAVE_FREE(x);
        RAVE_FREE(y);
        RAVE_FREE(uh);
        RAVE_FREE(vh);
        RAVE_FREE(vt1);
        RAVE_FREE(pointsTrigon);
        gsl_vector_free(uv);

        if(fitOk) return 1;
        else return 0;
}

max_stars_repo_name: adokter/radar

max_issues_repo_name: adokter/radar

max_forks_repo_name: adokter/radar


text:
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
//#include<ctype.h>
#include<string.h>
#include <gsl/gsl_fit.h>

#include "system_evol.h"
#include "initial_conditions.h"
#include "solver.h"
#include "chaotic_variables.h"

#define Distance(x, y) (sqrt((x) * (x) + (y) * (y)))
#define M 1 //masa de la particula
#define R sqrt(9) // define una region efectiva del potencial de radio R.
#define T 0
#define X 1
#define Y 2
#define VX 3
#define VY 4
#define En 5

extern double dt, t, phi; //declaradas en system_evol.c