text:
#pragma once

#include "SimpleEvaluator.h"
#include "BooleanDAG.h"
#include "Interface.h"
#include "SuggestionGenerator.h"
#include "OptimizationWrapper.h"
#include "MaxOptimizationWrapper.h"
#include "BoolNodeSimplifier.h"
#include "NumDebugger.h"
#include <map>
#include <set>

#ifndef _NOGSL
#include <gsl/gsl_vector.h>
#else
#include "FakeGSL.h"
#endif

#include "BoolBasedSampler.h"

class SuggestionGeneratorUsingMax: public SuggestionGenerator {
    Interface* interf;
    BooleanDAG* dag;
    map<string, int>& ctrls;
    FloatManager& fm;

    ActualEvaluators* actualEval;
    SmoothEvaluators* smoothEval;

    OptimizationWrapper* opt;
    MaxOptimizationWrapper* maxOpt;

    int minimizeNode;
    set<int> assertConstraints;

    vector<gsl_vector*> maxes;
    gsl_vector* tmp;

    int NUM_MAXES = 5;

    NumDebugger* debugger;
    const vector<vector<int>>& dependentInputs;

public:
    vector<gsl_vector*> dirGrads;
    SuggestionGeneratorUsingMax(BooleanDAG* _dag, FloatManager& _fm, Interface* _interface, map<string, int>& _ctrls, OptimizationWrapper* _opt, MaxOptimizationWrapper* _maxOpt, ActualEvaluators* _actualEval, SmoothEvaluators* _smoothEval, int ncontrols, const vector<vector<int>>& _dependentInputs, NumDebugger* _debugger): dag(_dag), fm(_fm), interf(_interface), ctrls(_ctrls), opt(_opt), maxOpt(_maxOpt), actualEval(_actualEval), smoothEval(_smoothEval), dependentInputs(_dependentInputs), debugger(_debugger) {

        minimizeNode = -1;

	    for (int i = 0; i < dag->size(); i++) { // TODO: this should also be set by caller class
	    	bool_node* n = (*dag)[i];
	    	if (n->type == bool_node::ASSERT && ((ASSERT_node*)n)->isHard()) {
	    		minimizeNode = i;
	    	} else if (n->type == bool_node::ASSERT || Util::isSqrt(n)) {
	    		assertConstraints.insert(i);
	    	} else if (n->type == bool_node::CTRL && n->getOtype() == OutType::BOOL) {
	    		assertConstraints.insert(i);
	    	}
	    }

        for (int i = 0; i < 20; i++) {
            dirGrads.push_back(gsl_vector_alloc(ncontrols));
        }

        for (int i = 0; i < NUM_MAXES; i++) {
            maxes.push_back(gsl_vector_alloc(ncontrols));
        }
        tmp = gsl_vector_alloc(ncontrols);
    }
    ~SuggestionGeneratorUsingMax() {
        for (int i = 0; i < dirGrads.size(); i++) {
            gsl_vector_free(dirGrads[i]);
        }

        for (int i = 0; i < maxes.size(); i++) {
            gsl_vector_free(maxes[i]);
        }
        gsl_vector_free(tmp);
    }

    pair<int, int> getSuggestion(const gsl_vector* state) {
        actualEval->run(state);

        vector<tuple<double, int, int>> s;

        for (auto it = interf->varsMapping.begin(); it != interf->varsMapping.end(); it++) {
            int nodeid = it->first;
            if (interf->hasValue(nodeid)) continue;

            bool_node* n = (*dag)[nodeid];
            bool hasArraccChild = Util::hasArraccChild(n);

            double dist = actualEval->dist(n->id);
            double cost = abs(dist);
            if (hasArraccChild && cost < 0.1) {
                return make_pair(nodeid, dist >= 0.0);
            }

            if (hasArraccChild) {
                cost = cost / 1000.0;
            }

            s.push_back(make_tuple(cost, nodeid, dist > 0));

        }
        cout << endl;

        sort(s.begin(), s.end());
        if (s.size() > 0) {
            return make_pair(get<1>(s[0]), get<2>(s[1]));
        } else {
            Assert(false, "weqhqp");
            return make_pair(-1, 0);
        }

    }

    pair<int, int> getUnsatSuggestion(const gsl_vector* state) {
        actualEval->run(state);
        double dassert;
        double d;
        vector<tuple<int, int>> conflict;
        for (auto it = assertConstraints.begin(); it != assertConstraints.end(); it++) {
            //dassert = actualEval->getErrorOnConstraint(*it);
            if (true) {
                // TODO: this is currently hacky
                const vector<int>& inputs = dependentInputs[*it];
                int minid = -1;
                double minerror  = 1e30;
                for (auto it1 = inputs.begin(); it1 != inputs.end(); it1++) {
                    if (interf->hasValue(*it1)) continue;
                    d = actualEval->dist(*it1);
                    //if (dassert < 0.0) {
                        cout << *it << " " << *it1 << " " << d << " " << Util::print(actualEval->grad(*it1)) << endl;
                    //}
                    if (d < minerror) {
                        minerror = d;
                        minid = *it1;
                    }
                }
                //Assert(minerror > 0.0, "something is wrong");
                if (minid != -1) {
                    cout << "Min " << *it << " " << minid << " " << minerror << endl;
                    conflict.push_back(make_tuple(*it, minid));
                }

                if ((*dag)[*it]->mother()->type == bool_node::NOT &&  minerror >= 0.2) {
                    break;
                }
            }

        }
        sort(conflict.begin(), conflict.end());
        reverse(conflict.begin(), conflict.end());
        Assert(conflict.size() > 0, "No failed asserts???");
        return make_pair(get<1>(conflict[0]), 1);

    }

    SClause* getUnsatClause(const gsl_vector* state, const gsl_vector* initState) {
        actualEval->run(state);
        double dassert;
        double d;
        vector<int> conflictingAsserts;
        for (auto it = assertConstraints.begin(); it != assertConstraints.end(); it++) {
            dassert = actualEval->getErrorOnConstraint(*it);
            if ( (*dag)[*it]->mother()->type != bool_node::NOT) {
                continue;
            }
            if (true) {
                if (dassert < 0.01) {
                    conflictingAsserts.push_back(*it);
                }
                if (dassert < -0.2) {
                    break;
                }
            }

        }

        SClause* clause = new SClause();

        actualEval->run(initState);
        for (auto it = conflictingAsserts.begin(); it != conflictingAsserts.end(); it++) {
                // TODO: this is currently hacky
                const vector<int>& inputs = dependentInputs[*it];
                int minid = -1;
                double minerror  = 1e30;
                for (auto it1 = inputs.begin(); it1 != inputs.end(); it1++) {
                    d = actualEval->dist(*it1);
                    cout << *it << " " << *it1 << " " << d << " " << Util::print(actualEval->grad(*it1)) << endl;
                    if (d < minerror) {
                        minerror = d;
                        minid = *it1;
                    }
                }
                //Assert(minerror > 0.0, "something is wrong");
                if (minid != -1 ) {
                    clause->add(*it, minid);
                }
        }

        return clause;

    }

    SClause* getUnsatClause(const gsl_vector* state) {
        actualEval->run(state);
        double dassert;
        double d;
        SClause* clause = new SClause();
        for (auto it = assertConstraints.begin(); it != assertConstraints.end(); it++) {
            dassert = actualEval->getErrorOnConstraint(*it);
            if (true) {
                // TODO: this is currently hacky
                const vector<int>& inputs = dependentInputs[*it];
                int minid = -1;
                double minerror  = 1e30;
                for (auto it1 = inputs.begin(); it1 != inputs.end(); it1++) {
                    d = actualEval->dist(*it1);
                    if (dassert < 0.0) {
                        cout << *it << " " << *it1 << " " << d << " " << Util::print(actualEval->grad(*it1)) << endl;
                    }
                    if (d < minerror) {
                        minerror = d;
                        minid = *it1;
                    }
                }
                //Assert(minerror > 0.0, "something is wrong");
                if (minid != -1 && (*dag)[*it]->mother()->type == bool_node::NOT &&  minerror >= -0.01) {
                    clause->add(*it, minid);
                }

                if ((*dag)[*it]->mother()->type == bool_node::NOT &&  minerror >= 0.2) {
                    break;
                }
            }

        }

        return clause;

    }

    BooleanDAG* getNewDag(BooleanDAG* origDag, int origNid, const gsl_vector* s) {
        BooleanDAG* newDag = origDag->clone();
        SimpleGradEvaluator* seval = actualEval->getEvaluator(origDag);
        seval->run(s);
        BoolNodeSimplifier simplifier(*newDag, fm);
        simplifier.process(*newDag, origNid, seval);
        DagOptim cse(*newDag, fm);
        cse.process(*newDag);
        newDag->cleanup();
        //newDag->lprint(cout);
        return newDag;
    }

    Predicate* make_pure(Predicate* p, const gsl_vector* s) {
        cout << "Making pure: " << p->print() << " at " << Util::print(s) << endl;
        if (p->isBasic()) {
            BasicPredicate* bp = (BasicPredicate*) p;
            BooleanDAG* newDag = getNewDag(bp->dag, bp->nid, s);
            // register newDag
            actualEval->addEvaluator(newDag);
            smoothEval->addEvaluator(newDag);
            int cid = -1;
            for (int i = newDag->size() - 1; i >= 0; i--) {
                bool_node* n = (*newDag)[i];
                if (n->type == bool_node::ASSERT) {
                    cid = n->mother()->id;
                    break;
                }
            }
            Assert(cid != -1, "No assert in new dag");
            return new BasicPredicate(newDag, cid, bp);
        } else {
            DiffPredicate* dp = (DiffPredicate*)p;
            Predicate* new_p1;
            Predicate* new_p2;
            if (dp->p1->isPure()) {
                new_p1 = dp->p1;
            } else {
                new_p1 = make_pure(dp->p1, s);
            }
            if (dp->p2->isPure()) {
                new_p2 = dp->p2;
            } else {
                new_p2 = make_pure(dp->p2, s);
            }
            int ncontrols = ctrls.size();
            if (ncontrols == 0) {
                ncontrols = 1;
            }
            return new DiffPredicate(new_p1, new_p2, ncontrols);
        }
    }

    int getDirectionsLevel0(const gsl_vector* state, double beta) {
        GradUtil::BETA = beta;
        GradUtil::ALPHA = -beta;
        smoothEval->run(state);

        vector<tuple<double, Predicate*, int>> s;

        for (auto it = interf->levelPredicates[0].begin(); it != interf->levelPredicates[0].end(); it++) {
            Predicate* p = *it;

            double dist = smoothEval->dist(p);
            double cost = abs(dist);
            if (cost >= 0.1) {
                s.push_back(make_tuple(cost, p, dist > 0));
            }
        }

        sort(s.begin(), s.end());
        cout << "Distances: " << endl;
        int end = 5;
        if (end > s.size()) {
            end = s.size();
        }
        for (int i = 0; i < end; i++) {
            cout << get<1>(s[i])->print() << " " << get<0>(s[i]) << " " << get<2>(s[i]) << endl;
        }

        int numDirs = 0;
        for (int i = 0; i < end; i++) {
            if (numDirs == 5) break;
            double dist = get<0>(s[i]);
            Predicate* p = get<1>(s[i]);
            int val = !get<2>(s[i]);
            if (dist < 0.1) continue;
            smoothEval->getErrorOnPredicate(p, val, dirGrads[numDirs]);
            cout << "G: " << Util::print(dirGrads[numDirs]) << endl;
            bool newDir = true;
            for (int j = 0; j < numDirs; j++) {
                if (Util::sameDir(dirGrads[j], dirGrads[numDirs])) {
                    newDir = false;
                    cout << "Not a new dir" << endl;
                    break;
                }
            }
            if (newDir) {
                numDirs++;
            }
        }
        Assert(numDirs > 0, "No dirs?");
        return numDirs;
    }

    void updateState(gsl_vector* state, int idx, double dx) {
        gsl_vector_set(state, idx, gsl_vector_get(state, idx) + dx);
    }

    void randomize(gsl_vector* tmp) {
        for (int i = 0; i < tmp->size; i++) {
            double r = -1.0 + (rand() % 20)/10.0;
            gsl_vector_set(tmp, i, r);
        }
    }

    int getMaxes(gsl_vector* state, double beta) {
        int numMaxes = 0;
        for (int i = 0; i < NUM_MAXES; i++) {
            randomize(tmp);
            gsl_vector_scale(tmp, 1.0/gsl_blas_dnrm2(tmp));
            gsl_vector_add(state, tmp);
            cout << Util::print(state) << endl;
            maxOpt->maximize(interf, state, tmp, assertConstraints, minimizeNode, beta, -1, i);
            gsl_vector_memcpy(maxes[i], maxOpt->getMinState());
            gsl_vector_sub(state, tmp);
            cout << Util::print(state) << endl;
            numMaxes++;
        }
        return numMaxes;
    }

    vector<tuple<double, Predicate*, int>> getDistancesLevel0(const gsl_vector* state) {
        actualEval->run(state);

        vector<tuple<double, Predicate*, int>> s;

        for (auto it = interf->levelPredicates[0].begin(); it != interf->levelPredicates[0].end(); it++) {
            Predicate* p = *it;

            double dist = actualEval->dist(p);
            double cost = abs(dist);
            s.push_back(make_tuple(cost, p, dist > 0));
        }

        sort(s.begin(), s.end());
        cout << "Distances: " << endl;
        int end = 5;
        if (end > s.size()) {
            end = s.size();
        }
        for (int i = 0; i < end; i++) {
            cout << get<1>(s[i])->print() << " " << get<0>(s[i]) << " " << get<2>(s[i]) << endl;
        }

        return s;
    }

    int chooseBeta(LocalState* ls) {
       double error1 = ls->errors[0];
       double error2 = ls->errors[1];
       double error3 = ls->errors[2];
       cout << "Errors: " << error1 << " " << error2 << " " << error3 << endl;

       if (error3 <= 0.0) {
            Assert(false, "Zero error?");
       }
       if (error3 < 1.0 || error2 <= 0.01) {
            return 2;
       }
       if (error2 < 5.0 || error1 <= 0.01) {
            return 1;
       } else {
            return 0;
       }
    }

    /*virtual IClause* getConflictClause(int level, LocalState* ls) {
        Assert(level == -1, "Other levels is not yet supported");
        int betaIdx = chooseBeta(ls);
        double beta;
        gsl_vector* state;
        if (betaIdx == 0) {
            beta = -1;
            state = ls->localSols[0];
        } else if (betaIdx == 1) {
            beta = -10;
            state = ls->localSols[1];
        } else if (betaIdx == 2) {
            beta = -50;
            state = ls->localSols[2];
        }
        cout << "Choosing beta = " << beta << endl;

        vector<tuple<double, Predicate*, int>> preds;
        actualEval->run(state);
        for (auto it = interf->levelPredicates[0].begin(); it != interf->levelPredicates[0].end(); it++) {
            Predicate* p = *it;
            cout << p->print() << " ";
            double dist = actualEval->dist(p);
            cout << dist << " ";
            cout << Util::print(actualEval->grad(p)) << endl;
            double cost = abs(dist);
            double norm = actualEval->grad_norm(p);
            if (norm < 0.01) continue;
            if (cost/norm < 0.01) continue;
            preds.push_back(make_tuple(cost/norm, p, dist >= 0));
        }
        sort(preds.begin(), preds.end());

        int numPreds = min(20, preds.size());
        IClause* c = new IClause(GradUtil::counter);
        for (int i = 0; i < numPreds; i++) {
            Predicate* p = get<1>(preds[i]);
            cout << p->print() << " ";
            if (!p->isPure()) {
                p = make_pure(p, state);
            }
            c->add(p, !get<2>(preds[i]));
        }
        cout << endl;

        if (c->size() == 0) {
            Assert(false, "Empty clause?");
        }
        return c;
    }*/

    virtual IClause* getConflictClause(int level, LocalState* ls) {
        int betaIdx = chooseBeta(ls);
        double beta;
        gsl_vector* state;
        if (betaIdx == 0) {
            beta = -1;
            state = ls->localSols[0];
        } else if (betaIdx == 1) {
            beta = -10;
            state = ls->localSols[1];
        } else if (betaIdx == 2) {
            beta = -50;
            state = ls->localSols[2];
        }
        cout << "Choosing beta = " << beta << endl;
        int numMaxes = getMaxes(state, beta);
        cout << "Num maxes: " << numMaxes << endl;

        map<int, vector<tuple<double, Predicate*, int>>> maxesToPredicates;
        map<Predicate*, vector<tuple<double, int, int>>> predicatesToMaxes;

        map<Predicate*, tuple<double, int>> lsPreds;

        actualEval->run(state);
        for (auto it = interf->levelPredicates[0].begin(); it != interf->levelPredicates[0].end(); it++) {
            Predicate* p = *it;
            double dist = actualEval->dist(p);
            double cost = abs(dist);
            double norm = actualEval->grad_norm(p);
            if (norm < 0.01) {
                cost = 1000.0;
            } else {
                cost = cost/norm;
            }
            lsPreds[p] = make_tuple(cost, dist >= 0);
        }

        for (int i = 0; i < numMaxes; i++) {
            actualEval->run(maxes[i]);
            vector<tuple<double, Predicate*, int>> preds;
            for (auto it = interf->levelPredicates[0].begin(); it != interf->levelPredicates[0].end(); it++) {
                Predicate* p = *it;
                double dist = actualEval->dist(p);
                double cost = abs(dist);
                double norm = actualEval->grad_norm(p);
                if (norm < 0.01) {
                    cost = 1000.0;
                } else {
                    cost = cost/norm;
                }
                auto& lsTup = lsPreds[p];
                bool same = false;
                if (get<1>(lsTup) == (dist > 0)) {
                    same = true;
                }
                if (get<0>(lsTup) <= 0.1) {
                    same = true;
                }
                if (cost <= 0.1) {
                    same = true;
                }
                if (!same) {
                    preds.push_back(make_tuple(cost, p, dist >= 0));
                    if (predicatesToMaxes.find(p) == predicatesToMaxes.end()) {
                        vector<tuple<double, int, int>> maxstates;
                        predicatesToMaxes[p] = maxstates;
                    }
                    predicatesToMaxes[p].push_back(make_tuple(cost, i, dist>=0));
                }

            }
            sort(preds.begin(), preds.end());
            maxesToPredicates[i] = (preds);
        }

        set<int> selected;

        IClause* c = new IClause(GradUtil::counter);

        for (int i = 0; i < numMaxes; i++) {
            if (selected.find(i) != selected.end()) continue;
            auto& l = maxesToPredicates[i];
            if (l.size() == 0) continue;
            //Assert(l.size() != 0, "No preds to eliminate max");
            auto& tup = l[0];
            Predicate* p = get<1>(tup);
            auto& m = predicatesToMaxes[p];
            selected.insert(i);
            for (int k = 0; k < m.size(); k++) {
                selected.insert(get<1>(m[k]));
            }
            if (!p->isPure()) {
                p = make_pure(p, state);
            }
            c->add(p, get<2>(tup));
        }

        if (c->size() == 0) {
            Assert(false, "Empty clause?");
        }

        // add predicates that satisfy all maxes and the current local state
        vector<tuple<double, Predicate*>> preds;
        actualEval->run(state);
        for (auto it = interf->levelPredicates[0].begin(); it != interf->levelPredicates[0].end(); it++) {
            Predicate* p = *it;
            auto& m = predicatesToMaxes[p];
            if (m.size() != 0) {
                continue;
            }
            double dist = actualEval->dist(p);
            double cost = abs(dist);
            double norm = actualEval->grad_norm(p);
            if (norm < 0.01) continue;
            preds.push_back(make_tuple(cost/norm, p));
        }
        sort(preds.begin(), preds.end());
        int numDirs = 0;
        vector<double> dists;
        dists.resize(dirGrads.size());
        for (int i = 0; i < preds.size(); i++) {
            if (numDirs >= dirGrads.size()) break;
            Predicate* p = get<1>(preds[i]);
            dists[numDirs] = actualEval->dist(p);
            gsl_vector_memcpy(dirGrads[numDirs], actualEval->grad(p));
            bool newDir = true;
            for (int j = 0; j < numDirs; j++) {
                if (Util::sameDir(dirGrads[j], dirGrads[numDirs])) {
                    if (dists[j] * dists[numDirs] >= 0.0) {
                        newDir = false;
                        //cout << "Not a new dir" << endl;
                        break;
                    }
                }
            }
            if (newDir) {
                cout << "New dir from ls: " << p->print() << " " << actualEval->dist(p) << " " << Util::print(actualEval->grad(p)) << endl;
                if (!p->isPure()) {
                    p = make_pure(p, state);
                }

                c->add(p, dists[numDirs] < 0);

                numDirs++;
            }
        }
        return c;

    }

    /*virtual IClause* getConflictClause(int level, LocalState* ls) {
        gsl_vector* state = ls->localSols[2];

        vector<tuple<double, Predicate*, int>> preds;
        actualEval->run(state);
        for (auto it = interf->levelPredicates[0].begin(); it != interf->levelPredicates[0].end(); it++) {
            Predicate* p = *it;
            double dist = actualEval->dist(p);
            double cost = abs(dist);
            double norm = actualEval->grad_norm(p);
            if (norm < 0.01) continue;
            if (cost/norm < 0.1) continue;
            preds.push_back(make_tuple(cost/norm, p, dist >= 0));
        }

        int numPreds = min(5, preds.size());
        IClause* c = new IClause(GradUtil::counter);
        for (int i = 0; i < numPreds; i++) {
            int ridx = rand() % preds.size();
            Predicate* p = get<1>(preds[ridx]);
            cout << p->print() << " " << actualEval->dist(p) << " " << Util::print(actualEval->grad(p)) << endl;
            if (!p->isPure()) {
                p = make_pure(p, state);
            }
            c->add(p, !get<2>(preds[ridx]));
        }

        if (c->size() == 0) {
            Assert(false, "Empty clause?");
        }
        return c;

    }*/

    /*virtual IClause* getConflictClause(int level, LocalState * ls) {
        Assert(level == -1, "Other levels is not yet supported");
        int betaIdx = chooseBeta(ls);
        double beta;
        gsl_vector* state;
        if (betaIdx == 0) {
            beta = -1;
            state = ls->localSols[0];
        } else if (betaIdx == 1) {
            beta = -10;
            state = ls->localSols[1];
        } else if (betaIdx == 2) {
            beta = -50;
            state = ls->localSols[2];
        }
        cout << "Choosing beta = " << beta << endl;
        int numMaxes = getMaxes(state, beta);
        cout << "Num maxes: " << numMaxes << endl;

        // collect close preds
        set<Predicate*> close_preds;

        vector<tuple<double, Predicate*>> preds;
        actualEval->run(state);
        for (auto it = interf->levelPredicates[0].begin(); it != interf->levelPredicates[0].end(); it++) {
            Predicate* p = *it;
            double dist = actualEval->dist(p);
            double cost = abs(dist);
            double norm = actualEval->grad_norm(p);
            if (norm < 0.01) continue;
            preds.push_back(make_tuple(cost/norm, p));
        }
        sort(preds.begin(), preds.end());
        int numDirs = 0;
        vector<double> dists;
        dists.resize(dirGrads.size());
        for (int i = 0; i < preds.size(); i++) {
            if (numDirs >= dirGrads.size()) break;
            Predicate* p = get<1>(preds[i]);
            dists[numDirs] = actualEval->dist(p);
            gsl_vector_memcpy(dirGrads[numDirs], actualEval->grad(p));
            bool newDir = true;
            for (int j = 0; j < numDirs; j++) {
                if (Util::sameDir(dirGrads[j], dirGrads[numDirs])) {
                    if (dists[j] * dists[numDirs] >= 0.0) {
                        newDir = false;
                        //cout << "Not a new dir" << endl;
                        break;
                    }
                }
            }
            if (newDir) {
                cout << "New dir from ls: " << p->print() << " " << actualEval->dist(p) << " " << Util::print(actualEval->grad(p)) << endl;
                numDirs++;
                if (p->isPure()) {
                    close_preds.insert(p);
                } else {
                    p = make_pure(p, state);
                    close_preds.insert(p);
                }
            }
        }

        for (int k = 0; k < numMaxes; k++) {
            preds.clear();
            actualEval->run(maxes[k]);
            for (auto it = interf->levelPredicates[0].begin(); it != interf->levelPredicates[0].end(); it++) {
                Predicate* p = *it;

                double dist = actualEval->dist(p);
                double cost = abs(dist);
                double norm = actualEval->grad_norm(p);
                if (norm < 0.01) continue;
                preds.push_back(make_tuple(cost/norm, p));
            }
            sort(preds.begin(), preds.end());
            int numDirs = 0;
            for (int i = 0; i < preds.size(); i++) {
                if (numDirs >= 20) break;
                Predicate* p = get<1>(preds[i]);
                dists[numDirs] = actualEval->dist(p);
                gsl_vector_memcpy(dirGrads[numDirs], actualEval->grad(p));
                bool newDir = true;
                for (int j = 0; j < numDirs; j++) {
                    if (Util::sameDir(dirGrads[j], dirGrads[numDirs])) {
                        if (dists[j] * dists[numDirs] >= 0.0) {
                            newDir = false;
                            //cout << "Not a new dir" << endl;
                            break;
                        }
                    }
                }
                if (newDir) {
                    cout << "New dir from max: " << p->print() << " " << actualEval->dist(p) << " " << Util::print(actualEval->grad(p)) << endl;
                    numDirs++;
                    if (p->isPure()) {
                        close_preds.insert(p);
                    } else {
                        p = make_pure(p, maxes[k]);
                        close_preds.insert(p);
                    }
                }
            }
        }
        // collect common conditions satisfied by state, maxes
        map<Predicate*, vector<tuple<double, int>>> common_preds;

        actualEval->run(state);
        for (auto it = close_preds.begin(); it != close_preds.end(); it++) {
            Predicate* p = *it;
            cout << "P: " << p->print() << " ";
            double dist = actualEval->dist(p);
            cout << dist << " ";
            cout << Util::print(actualEval->grad(p)) << endl;
            double cost = abs(dist);
            int val = dist >= 0;
            vector<tuple<double, int>> costValPairs;
            costValPairs.push_back(make_tuple(cost, val));
            common_preds[p] = costValPairs;
        }

        for (int i = 0; i < numMaxes; i++) {
            actualEval->run(maxes[i]);
            for (auto it = close_preds.begin(); it != close_preds.end(); it++) {
                Predicate* p = *it;

                double dist = actualEval->dist(p);
                double cost = abs(dist);
                int val = dist >= 0;
                common_preds[p].push_back(make_tuple(cost, val));
            }
        }

        vector<tuple<double, Predicate*, int>> ss;

        for (auto it = common_preds.begin(); it != common_preds.end(); it++) {
            Predicate* p = it->first;

            vector<tuple<double, int>>& pairs = common_preds[p];
            int val = -1;
            double cost = 1e10;

            cout << "Predicate " << p->print() << " " ;

            for (int i = 0; i < pairs.size(); i++) {
                auto& tup = pairs[i];
                double tcost = get<0>(tup);
                double tval = get<1>(tup);
                cout << "(" << tval << "," << tcost << ")" << " ";
                if (tcost < cost) {
                    cost = tcost;
                }
                if (tcost > 0.01) {
                    if (val == -1) {
                        val = tval;
                    } else if (val != tval) {
                        val = -2;
                        break;
                    }
                }
            }
            cout << "Final: " << val << "," << cost << endl;
            if (val == 0 || val == 1) {
                ss.push_back(make_tuple(cost, p, val));
            }
        }
        sort(ss.begin(), ss.end());

        cout << "Num common preds: " << ss.size() << endl;

        int numPreds = min(20, ss.size());
        IClause* c = new IClause(GradUtil::counter);
        for (int i = 0; i < numPreds; i++) {
            c->add(get<1>(ss[i]), !get<2>(ss[i]));
        }

        if (c->size() == 0) {
            Assert(false, "Empty clause?");
        }
        return c;
    }*/

};

max_stars_repo_name: natebragg/sketch-backend

max_issues_repo_name: natebragg/sketch-backend

max_forks_repo_name: natebragg/sketch-backend


text:
#include <math.h>
#include <string.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

extern void getSeries(int dataPoints, double dt, int nComps, int nTypes, int* epsilon, int capacity, int* curState, int rng_seed, int* output);

ext: h

max_stars_repo_name: akononovicius/compartmental-voter-model

max_issues_repo_name: akononovicius/compartmental-voter-model

max_forks_repo_name: akononovicius/compartmental-voter-model


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

/* Factorise a symmetric matrix A into
 *
 * A = Q T Q'
 *
 * where Q is orthogonal and T is symmetric tridiagonal.  Only the
 * diagonal and lower triangular part of A is referenced and modified.
 *
 * On exit, T is stored in the diagonal and first subdiagonal of
 * A. Since T is symmetric the upper diagonal is not stored.
 *
 * Q is stored as a packed set of Householder transformations in the
 * lower triangular part of the input matrix below the first subdiagonal.
 *
 * The full matrix for Q can be obtained as the product
 *
 *       Q = Q_1 Q_2 ... Q_(N-2)
 *
 * where
 *
 *       Q_i = (I - tau_i * v_i * v_i')
 *
 * and where v_i is a Householder vector
 *
 *       v_i = [0, ... , 0, 1, A(i+1,i), A(i+2,i), ... , A(N,i)]
 *
 * This storage scheme is the same as in LAPACK.  See LAPACK's
 * ssytd2.f for details.
 *
 * See Golub & Van Loan, "Matrix Computations" (3rd ed), Section 8.3
 *
 * Note: this description uses 1-based indices. The code below uses
 * 0-based indices
 */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>

#include <gsl/gsl_linalg.h>

int
gsl_linalg_symmtd_decomp (gsl_matrix * A, gsl_vector * tau)
{
  if (A->size1 != A->size2)
    {
      GSL_ERROR ("symmetric tridiagonal decomposition requires square matrix",
                 GSL_ENOTSQR);
    }
  else if (tau->size + 1 != A->size1)
    {
      GSL_ERROR ("size of tau must be N-1", GSL_EBADLEN);
    }
  else
    {
      const size_t N = A->size1;
      size_t i;

      for (i = 0 ; i < N - 2; i++)
        {
          gsl_vector_view v = gsl_matrix_subcolumn (A, i, i + 1, N - i - 1);
          double tau_i = gsl_linalg_householder_transform (&v.vector);

          /* Apply the transformation H^T A H to the remaining columns */

          if (tau_i != 0.0)
            {
              gsl_matrix_view m = gsl_matrix_submatrix (A, i + 1, i + 1, N - i - 1, N - i - 1);
              double ei = gsl_vector_get(&v.vector, 0);
              gsl_vector_view x = gsl_vector_subvector (tau, i, N - i - 1);

              gsl_vector_set (&v.vector, 0, 1.0);

              /* x = tau * A * v */
              gsl_blas_dsymv (CblasLower, tau_i, &m.matrix, &v.vector, 0.0, &x.vector);

              /* w = x - (1/2) tau * (x' * v) * v  */
              {
                double xv, alpha;
                gsl_blas_ddot(&x.vector, &v.vector, &xv);
                alpha = -0.5 * tau_i * xv;
                gsl_blas_daxpy(alpha, &v.vector, &x.vector);
              }

              /* apply the transformation A = A - v w' - w v' */
              gsl_blas_dsyr2(CblasLower, -1.0, &v.vector, &x.vector, &m.matrix);

              gsl_vector_set (&v.vector, 0, ei);
            }

          gsl_vector_set (tau, i, tau_i);
        }

      return GSL_SUCCESS;
    }
}

/*  Form the orthogonal matrix Q from the packed QR matrix */

int
gsl_linalg_symmtd_unpack (const gsl_matrix * A,
                          const gsl_vector * tau,
                          gsl_matrix * Q,
                          gsl_vector * diag,
                          gsl_vector * sdiag)
{
  if (A->size1 !=  A->size2)
    {
      GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
    }
  else if (tau->size + 1 != A->size1)
    {
      GSL_ERROR ("size of tau must be (matrix size - 1)", GSL_EBADLEN);
    }
  else if (Q->size1 != A->size1 || Q->size2 != A->size1)
    {
      GSL_ERROR ("size of Q must match size of A", GSL_EBADLEN);
    }
  else if (diag->size != A->size1)
    {
      GSL_ERROR ("size of diagonal must match size of A", GSL_EBADLEN);
    }
  else if (sdiag->size + 1 != A->size1)
    {
      GSL_ERROR ("size of subdiagonal must be (matrix size - 1)", GSL_EBADLEN);
    }
  else
    {
      const size_t N = A->size1;
      gsl_vector_const_view d = gsl_matrix_const_diagonal(A);;
      gsl_vector_const_view sd = gsl_matrix_const_subdiagonal(A, 1);;
      size_t i;

      /* Initialize Q to the identity */

      gsl_matrix_set_identity (Q);

      for (i = N - 2; i-- > 0;)
        {
          gsl_vector_const_view h = gsl_matrix_const_subcolumn (A, i, i + 1, N - i - 1);
          double ti = gsl_vector_get (tau, i);
          gsl_matrix_view m = gsl_matrix_submatrix (Q, i + 1, i + 1, N - i - 1, N - i - 1);
          gsl_vector_view work = gsl_vector_subvector(diag, 0, N - i - 1);

          gsl_linalg_householder_left (ti, &h.vector, &m.matrix, &work.vector);
        }

      /* copy diagonal into diag */
      gsl_vector_memcpy(diag, &d.vector);

      /* copy subdiagonal into sd */
      gsl_vector_memcpy(sdiag, &sd.vector);

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_symmtd_unpack_T (const gsl_matrix * A,
                            gsl_vector * diag,
                            gsl_vector * sdiag)
{
  if (A->size1 !=  A->size2)
    {
      GSL_ERROR ("matrix A must be square", GSL_ENOTSQR);
    }
  else if (diag->size != A->size1)
    {
      GSL_ERROR ("size of diagonal must match size of A", GSL_EBADLEN);
    }
  else if (sdiag->size + 1 != A->size1)
    {
      GSL_ERROR ("size of subdiagonal must be (matrix size - 1)", GSL_EBADLEN);
    }
  else
    {
      const size_t N = A->size1;
      gsl_vector_const_view d = gsl_matrix_const_diagonal(A);;
      gsl_vector_const_view sd = gsl_matrix_const_subdiagonal(A, 1);;

      /* copy diagonal into diag */
      gsl_vector_memcpy(diag, &d.vector);

      /* copy subdiagonal into sd */
      gsl_vector_memcpy(sdiag, &sd.vector);

      return GSL_SUCCESS;
    }
}

ext: c

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
/**
 * This file is part of the "libterminal" project
 *   Copyright (c) 2019-2020 Christian Parpart <christian@parpart.family>
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
#pragma once

#include <terminal/Color.h>
#include <terminal/RenderBuffer.h>
#include <terminal/Screen.h>

#include <terminal_renderer/BoxDrawingRenderer.h>
#include <terminal_renderer/FontDescriptions.h>
#include <terminal_renderer/RenderTarget.h>
#include <terminal_renderer/TextureAtlas.h>

#include <text_shaper/font.h>
#include <text_shaper/shaper.h>

#include <crispy/FNV.h>
#include <crispy/LRUCache.h>
#include <crispy/point.h>
#include <crispy/size.h>

#include <unicode/convert.h>
#include <unicode/run_segmenter.h>

#include <gsl/span>
#include <gsl/span_ext>

#include <functional>
#include <list>
#include <memory>
#include <unordered_map>
#include <vector>

namespace terminal::renderer
{

std::unique_ptr<text::font_locator> createFontLocator(FontLocatorEngine _engine);

struct FontKeys
{
    text::font_key regular;
    text::font_key bold;
    text::font_key italic;
    text::font_key boldItalic;
    text::font_key emoji;
};

/// Text Rendering Pipeline
class TextRenderer: public Renderable
{
  public:
    TextRenderer(GridMetrics const& gridMetrics,
                 text::shaper& textShaper,
                 FontDescriptions& fontDescriptions,
                 FontKeys const& fontKeys);

    void setRenderTarget(RenderTarget& renderTarget, DirectMappingAllocator& directMappingAllocator) override;
    void setTextureAtlas(TextureAtlas& atlas) override;

    void inspect(std::ostream& _textOutput) const override;

    void clearCache() override;

    void updateFontMetrics();

    void setPressure(bool _pressure) noexcept { pressure_ = _pressure; }

    /// Must be invoked before a new terminal frame is rendered.
    void beginFrame();

    /// Renders a given terminal's grid cell that has been
    /// transformed into a RenderCell.
    void renderCell(RenderCell const& _cell);

    /// Must be invoked when rendering the terminal's text has finished for this frame.
    void endFrame();

  private:
    void initializeDirectMapping();

    /// Puts a sequence of codepoints that belong to the same grid cell at @p _pos
    /// at the end of the currently filled line.
    void appendCellTextToClusterGroup(std::u32string const& _codepoints, TextStyle _style, RGBColor _color);

    /// Gets the text shaping result of the current text cluster group
    text::shape_result const& getOrCreateCachedGlyphPositions(crispy::StrongHash hash);
    text::shape_result createTextShapedGlyphPositions();
    text::shape_result shapeTextRun(unicode::run_segmenter::range const& _run);
    void flushTextClusterGroup();

    AtlasTileAttributes const* getOrCreateRasterizedMetadata(crispy::StrongHash const& hash,
                                                             text::glyph_key const& glyphKey,
                                                             unicode::PresentationStyle presentationStyle);

    /**
     * Creates (and rasterizes) a single glyph and returns its
     * render tile attributes required for the render step.
     */
    std::optional<TextureAtlas::TileCreateData> createSlicedRasterizedGlyph(
        atlas::TileLocation tileLocation,
        text::glyph_key const& id,
        unicode::PresentationStyle presentation,
        crispy::StrongHash const& hash);

    std::optional<TextureAtlas::TileCreateData> createRasterizedGlyph(
        atlas::TileLocation tileLocation, text::glyph_key const& id, unicode::PresentationStyle presentation);

    crispy::Point applyGlyphPositionToPen(crispy::Point pen,
                                          AtlasTileAttributes const& tileAttributes,
                                          text::glyph_position const& gpos) const noexcept;

    void renderRasterizedGlyph(crispy::Point targetSurfacePosition,
                               RGBAColor glyphColor,
                               AtlasTileAttributes const& attributes);

    // general properties
    //
    FontDescriptions& fontDescriptions_;
    FontKeys const& fonts_;

    // performance optimizations
    //
    bool pressure_ = false;

    using ShapingResultCache = crispy::StrongLRUHashtable<text::shape_result>;
    using ShapingResultCachePtr = ShapingResultCache::Ptr;

    ShapingResultCachePtr textShapingCache_;
    // TODO: make unique_ptr, get owned, export cref for other users in Renderer impl.
    text::shaper& textShaper_;

    DirectMapping _directMapping {};

    // Maps from glyph index to tile index.
    std::vector<uint32_t> _directMappedGlyphKeyToTileIndex {};

    bool isGlyphDirectMapped(text::glyph_key const& glyph) const noexcept
    {
        return _directMapping                  // Is direct mapping enabled?
               && glyph.font == fonts_.regular // Only regular font is direct-mapped for now.
               && glyph.index.value < _directMappedGlyphKeyToTileIndex.size()
               && _directMappedGlyphKeyToTileIndex[glyph.index.value] != 0;
    }

    // sub-renderer
    //
    BoxDrawingRenderer boxDrawingRenderer_;

    // work-data for the current text cluster group
    struct TextClusterGroup
    {
        // pen-start position of this text group
        crispy::Point initialPenPosition {};

        // uniform text style for this text group
        TextStyle style = TextStyle::Invalid;

        // uniform text color for this text group
        RGBColor color {};

        // codepoints within this text group with
        // uniform unicode properties (script, language, direction).
        std::vector<char32_t> codepoints;

        // cluster indices for each codepoint
        std::vector<unsigned> clusters;

        // number of grid cells processed
        int cellCount = 0; // FIXME: EA width vs actual cells

        void resetAndMovePenForward(int penIncrementInX)
        {
            codepoints.clear();
            clusters.clear();
            cellCount = 0;
            initialPenPosition.x += penIncrementInX;
        }
    };
    TextClusterGroup textClusterGroup_ {};

    bool textStartFound_ = false;
    bool forceCellGroupSplit_ = false;
};

} // namespace terminal::renderer

max_stars_repo_name: dandycheung/contour

max_issues_repo_name: dandycheung/contour

max_forks_repo_name: dandycheung/contour


text:
/* randist/gamma.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 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.
 */

#include <config.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

static double gamma_large (const gsl_rng * r, const double a);
static double gamma_frac (const gsl_rng * r, const double a);

/* The Gamma distribution of order a>0 is defined by:

   p(x) dx = {1 / \Gamma(a) b^a } x^{a-1} e^{-x/b} dx

   for x>0.  If X and Y are independent gamma-distributed random
   variables of order a1 and a2 with the same scale parameter b, then
   X+Y has gamma distribution of order a1+a2.

   The algorithms below are from Knuth, vol 2, 2nd ed, p. 129. */

double
gsl_ran_gamma_knuth (const gsl_rng * r, const double a, const double b)
{
  /* assume a > 0 */
  unsigned int na = floor (a);

  if (a == na)
    {
      return b * gsl_ran_gamma_int (r, na);
    }
  else if (na == 0)
    {
      return b * gamma_frac (r, a);
    }
  else
    {
      return b * (gsl_ran_gamma_int (r, na) + gamma_frac (r, a - na)) ;
    }
}

double
gsl_ran_gamma_int (const gsl_rng * r, const unsigned int a)
{
  if (a < 12)
    {
      unsigned int i;
      double prod = 1;

      for (i = 0; i < a; i++)
        {
          prod *= gsl_rng_uniform_pos (r);
        }

      /* Note: for 12 iterations we are safe against underflow, since
         the smallest positive random number is O(2^-32). This means
         the smallest possible product is 2^(-12*32) = 10^-116 which
         is within the range of double precision. */

      return -log (prod);
    }
  else
    {
      return gamma_large (r, (double) a);
    }
}

static double
gamma_large (const gsl_rng * r, const double a)
{
  /* Works only if a > 1, and is most efficient if a is large

     This algorithm, reported in Knuth, is attributed to Ahrens.  A
     faster one, we are told, can be found in: J. H. Ahrens and
     U. Dieter, Computing 12 (1974) 223-246.  */

  double sqa, x, y, v;
  sqa = sqrt (2 * a - 1);
  do
    {
      do
        {
          y = tan (M_PI * gsl_rng_uniform (r));
          x = sqa * y + a - 1;
        }
      while (x <= 0);
      v = gsl_rng_uniform (r);
    }
  while (v > (1 + y * y) * exp ((a - 1) * log (x / (a - 1)) - sqa * y));

  return x;
}

static double
gamma_frac (const gsl_rng * r, const double a)
{
  /* This is exercise 16 from Knuth; see page 135, and the solution is
     on page 551.  */

  double p, q, x, u, v;
  p = M_E / (a + M_E);
  do
    {
      u = gsl_rng_uniform (r);
      v = gsl_rng_uniform_pos (r);

      if (u < p)
        {
          x = exp ((1 / a) * log (v));
          q = exp (-x);
        }
      else
        {
          x = 1 - log (v);
          q = exp ((a - 1) * log (x));
        }
    }
  while (gsl_rng_uniform (r) >= q);

  return x;
}

double
gsl_ran_gamma_pdf (const double x, const double a, const double b)
{
  if (x < 0)
    {
      return 0 ;
    }
  else if (x == 0)
    {
      if (a == 1)
        return 1/b ;
      else
        return 0 ;
    }
  else if (a == 1)
    {
      return exp(-x/b)/b ;
    }
  else
    {
      double p;
      double lngamma = gsl_sf_lngamma (a);
      p = exp ((a - 1) * log (x/b) - x/b - lngamma)/b;
      return p;
    }
}

/* New version based on Marsaglia and Tsang, "A Simple Method for
 * generating gamma variables", ACM Transactions on Mathematical
 * Software, Vol 26, No 3 (2000), p363-372.
 *
 * Implemented by J.D.Lamb@btinternet.com, minor modifications for GSL
 * by Brian Gough
 */

double
gsl_ran_gamma_mt (const gsl_rng * r, const double a, const double b)
{
  return gsl_ran_gamma (r, a, b);
}

double
gsl_ran_gamma (const gsl_rng * r, const double a, const double b)
{
  /* assume a > 0 */

  if (a < 1)
    {
      double u = gsl_rng_uniform_pos (r);
      return gsl_ran_gamma (r, 1.0 + a, b) * pow (u, 1.0 / a);
    }

  {
    double x, v, u;
    double d = a - 1.0 / 3.0;
    double c = (1.0 / 3.0) / sqrt (d);

    while (1)
      {
        do
          {
            x = gsl_ran_gaussian_ziggurat (r, 1.0);
            v = 1.0 + c * x;
          }
        while (v <= 0);

        v = v * v * v;
        u = gsl_rng_uniform_pos (r);

        if (u < 1 - 0.0331 * x * x * x * x)
          break;

        if (log (u) < 0.5 * x * x + d * (1 - v + log (v)))
          break;
      }

    return b * d * v;
  }
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


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

#include <config.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_multiset.h>

size_t
gsl_multiset_n (const gsl_multiset * c)
{
  return c->n ;
}

size_t
gsl_multiset_k (const gsl_multiset * c)
{
  return c->k ;
}

size_t *
gsl_multiset_data (const gsl_multiset * c)
{
  return c->data ;
}

int
gsl_multiset_valid (gsl_multiset * c)
{
  const size_t n = c->n ;
  const size_t k = c->k ;

  size_t i, j ;

  for (i = 0; i < k; i++)
    {
      const size_t ci = c->data[i];

      if (ci >= n)
        {
          GSL_ERROR("multiset index outside range", GSL_FAILURE) ;
        }

      for (j = 0; j < i; j++)
        {
          if (c->data[j] > ci)
            {
              GSL_ERROR("multiset indices not in increasing order",
                        GSL_FAILURE) ;
            }
        }
    }

  return GSL_SUCCESS;
}

int
gsl_multiset_next (gsl_multiset * c)
{
  /* Replaces c with the next multiset (in the standard lexicographical
   * ordering).  Returns GSL_FAILURE if there is no next multiset.
   */
  const size_t n = c->n;
  const size_t k = c->k;
  size_t *data = c->data;
  size_t i;

  if(k == 0)
    {
      return GSL_FAILURE;
    }
  i = k - 1;

  while(i > 0 && data[i] == n-1)
    {
      --i;
    }

  if (i == 0 && data[0] == n-1)
    {
      return GSL_FAILURE;
    }

  ++data[i];

  while(i < k-1)
    {
      data[i+1] = data[i];
      ++i;
    }

  return GSL_SUCCESS;
}

int
gsl_multiset_prev (gsl_multiset * c)
{
  /* Replaces c with the previous multiset (in the standard
   * lexicographical ordering).  Returns GSL_FAILURE if there is no
   * previous multiset.
   */
  const size_t n = c->n;
  const size_t k = c->k;
  size_t *data = c->data;
  size_t i;

  if(k == 0)
    {
      return GSL_FAILURE;
    }
  i = k - 1;

  while(i > 0 && data[i-1] == data[i])
    {
      --i;
    }

  if(i == 0 && data[i] == 0)
    {
      return GSL_FAILURE;
    }

  data[i]--;

  if (data[i] < n-1)
    {
      while (i < k-1) {
        data[++i] = n - 1;
      }
    }

  return GSL_SUCCESS;
}

int
gsl_multiset_memcpy (gsl_multiset * dest, const gsl_multiset * src)
{
   const size_t src_n = src->n;
   const size_t src_k = src->k;
   const size_t dest_n = dest->n;
   const size_t dest_k = dest->k;

   if (src_n != dest_n || src_k != dest_k)
     {
       GSL_ERROR ("multiset lengths are not equal", GSL_EBADLEN);
     }

   {
     size_t j;

     for (j = 0; j < src_k; j++)
       {
         dest->data[j] = src->data[j];
       }
   }

   return GSL_SUCCESS;
}

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/* Siconos is a program dedicated to modeling, simulation and control
 * of non smooth dynamical systems.
 *
 * Copyright 2021 INRIA.
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 * http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
*/

#ifndef SiconosLAPACKE_H
#define SiconosLAPACKE_H

// IWYU pragma: private, include "SiconosLapack.h"
//#include "SiconosBlas.h"
#include <stdlib.h>
#include <assert.h>
#include <stdio.h>

// -------- Headers and routines naming conventions for the different Lapack implementations --------

// --- Intel MKL Header ---
#if defined(HAS_MKL_LAPACKE)
#include <mkl_lapacke.h>
#else
// Standard lapacke header
#include <lapacke.h>
#endif

// Name of the routines
#define LAPACK_NAME(N) LAPACKE_##N

#define LA_TRANS 'T'
#define LA_NOTRANS 'N'
#define LA_UP 'U'
#define LA_LO 'L'
#define LA_NONUNIT 'N'
#define LA_UNIT 'U'

#define INTEGER(X) X
#define INTEGERP(X) X
#define CHAR(X) X

#ifndef lapack_int
#define lapack_int int
#endif

// --- DGESVD ---
#if defined(HAS_LAPACK_dgesvd)
#define WRAP_DGESVD(F,A1,A2,A3,A4,A5,A6,A7,A8,A9,A10,A11,A12,INFO)      \
  INFO = F(CblasColMajor,A1,A2,A3,A4,A5,A6,A7,A8,A9,A10,A11,A12)
#else
#define WRAP_DGESVD(F,A1,A2,A3,A4,A5,A6,A7,A8,A9,A10,A11,A12,INFO)      \
  fprintf(stderr, "Your lapack version misses dgesvd function.\n");
#endif

// --- DGETRS ---
#define WRAP_DGETRS(F,A1,A2,A3,A4,A5,A6,A7,A8,INFO)   \
  INFO = F(CblasColMajor,A1,A2,A3,A4,A5,A6,A7,A8)

// --- DPOTRS ---
#define WRAP_DPOTRS(F,A1,A2,A3,A4,A5,A6,A7,INFO) \
  INFO = F(CblasColMajor,A1,A2,A3,A4,A5,A6,A7)

// --- DSYTRS ---
#define WRAP_DSYTRS(F,A1,A2,A3,A4,A5,A6,A7,A8,INFO) \
  INFO = F(CblasColMajor,A1,A2,A3,A4,A5,A6,A7,A8)

// --- DGESV ---
#define WRAP_DGESV(F,A1,A2,A3,A4,A5,A6,A7,INFO)   \
  INFO = F(CblasColMajor,A1,A2,A3,A4,A5,A6,A7)

// --- DPOSV ---
#define WRAP_DPOSV(F,A1,A2,A3,A4,A5,A6,A7,INFO) \
  INFO = F(CblasColMajor,A1,A2,A3,A4,A5,A6,A7)

// --- DGELS ---
#if defined(HAS_LAPACK_dgels)
#define WRAP_DGELS(F,A1,A2,A3,A4,A5,A6,A7,A8,INFO)    \
  INFO = F(CblasColMajor,A1,A2,A3,A4,A5,A6,A7,A8)
#else
#define WRAP_DGELS(F,A1,A2,A3,A4,A5,A6,A7,A8,INFO)                      \
  fprintf(stderr, "Your lapack version misses dgels function.\n");
#endif

// --- DGETRI ---
#define WRAP_DGETRI(F,A1,A2,A3,A4,INFO)         \
  INFO = F(CblasColMajor,A1,A2,A3,A4)

// --- DGETRF ---
#define WRAP_DGETRF(F,A1,A2,A3,A4,A5,INFO)  \
  INFO = F(CblasColMajor,A1,A2,A3,A4,A5)

// --- DPOTRF ---
#define WRAP_DPOTRF(F,A1,A2,A3,A4,INFO)  \
  INFO = F(CblasColMajor,A1,A2,A3,A4)

// --- DSYTRF ---
#define WRAP_DSYTRF(F,A1,A2,A3,A4,A5,INFO)      \
  INFO = F(CblasColMajor,A1,A2,A3,A4,A5)

// --- DTRTRS ---
#if defined(HAS_LAPACK_dtrtrs)
#define WRAP_DTRTRS(F,A1,A2,A3,A4,A5,A6,A7,A8,A9,INFO)  \
  INFO = F(CblasColMajor,A1,A2,A3,A4,A5,A6,A7,A8,A9)
#else
#define WRAP_DTRTRS(F,A1,A2,A3,A4,A5,A6,A7,A8,A9,INFO)                \
  fprintf(stderr, "Your lapack version misses dtrtrs function.\n");
#endif

#endif // SICONOSLAPACKE_H

max_stars_repo_name: BuildJet/siconos

max_issues_repo_name: BuildJet/siconos

max_forks_repo_name: BuildJet/siconos


text:
/*
* Copyright (c) 2006 Rudi Cilibrasi, Rulers of the RHouse
* All rights reserved.     cilibrar@cilibrar.com
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
*     * Redistributions of source code must retain the above copyright
*       notice, this list of conditions and the following disclaimer.
*     * 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.
*     * Neither the name of the RHouse nor the
*       names of its contributors may be used to endorse or promote products
*       derived from this software without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE RULERS AND CONTRIBUTORS "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 THE RULERS AND 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.
*/
#ifndef __CLOUTPUT_H
#define __CLOUTPUT_H

#include <math.h>

#include <gsl/gsl_blas.h>

/*! \file cloutput.h */

gsl_matrix *complearn_svd_project(gsl_matrix *a);
GString *complearn_matrix_prettyprint_text(LabeledMatrix *inp);
GString *complearn_matrix_prettyprint_nex(LabeledMatrix *inp, GString *treeblock);

struct CLDistMatrix {
  char *fileverstr;
  char *cllibver;
  char *username;
  char *hostname;
  char *title;
  char *compressor;
  char *creationTime; // Seconds since the epoch
  char * *cmds;
  char * *cmdtimes; // Seconds since the epoch
  LabeledMatrix *m;
  GArray *cmdsa;
  GArray *cmdtimesa;
  GArray *numa;
  GArray *labels1a;
  GArray *labels2a;
};

struct CLDistMatrix *complearn_read_clb_dist_matrix(const GString *db);
GArray *complearn_get_nexus_labels(const GString *db);
gsl_matrix *complearn_get_nexus_distance_matrix(const GString *db);
GString *complearn_matrix_prettyprint_clb(LabeledMatrix *result);
LabeledMatrix *complearn_load_nexus_matrix(const GString *inp);
LabeledMatrix *complearn_load_any_matrix(const GString *inp);
gboolean complearn_is_nexus_file(const GString *db);
GString *complearn_get_nexus_tree_block(const GString *db);
gboolean complearn_is_text_matrix(const GString *db);
LabeledMatrix *complearn_load_text_matrix(const GString *inp);

#endif

max_stars_repo_name: rudi-cilibrasi/classic-complearn

max_issues_repo_name: rudi-cilibrasi/classic-complearn

max_forks_repo_name: rudi-cilibrasi/classic-complearn


text:
/* specfunc/gamma_inc.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_erf.h>
#include <gsl/gsl_sf_exp.h>
#include <gsl/gsl_sf_log.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_expint.h>

#include "error.h"

/* The dominant part,
 * D(a,x) := x^a e^(-x) / Gamma(a+1)
 */
static
int
gamma_inc_D(const double a, const double x, gsl_sf_result * result)
{
  if(a < 10.0) {
    double lnr;
    gsl_sf_result lg;
    gsl_sf_lngamma_e(a+1.0, &lg);
    lnr = a * log(x) - x - lg.val;
    result->val = exp(lnr);
    result->err = 2.0 * GSL_DBL_EPSILON * (fabs(lnr) + 1.0) * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    gsl_sf_result gstar;
    gsl_sf_result ln_term;
    double term1;
    if (x < 0.5*a) {
      double u = x/a;
      double ln_u = log(u);
      ln_term.val = ln_u - u + 1.0;
      ln_term.err = (fabs(ln_u) + fabs(u) + 1.0) * GSL_DBL_EPSILON;
    } else {
      double mu = (x-a)/a;
      gsl_sf_log_1plusx_mx_e(mu, &ln_term);  /* log(1+mu) - mu */
    };
    gsl_sf_gammastar_e(a, &gstar);
    term1 = exp(a*ln_term.val)/sqrt(2.0*M_PI*a);
    result->val  = term1/gstar.val;
    result->err  = 2.0 * GSL_DBL_EPSILON * (fabs(a*ln_term.val) + 1.0) * fabs(result->val);
    result->err += gstar.err/fabs(gstar.val) * fabs(result->val);
    return GSL_SUCCESS;
  }

}

/* P series representation.
 */
static
int
gamma_inc_P_series(const double a, const double x, gsl_sf_result * result)
{
  const int nmax = 5000;

  gsl_sf_result D;
  int stat_D = gamma_inc_D(a, x, &D);

  double sum  = 1.0;
  double term = 1.0;
  int n;
  for(n=1; n<nmax; n++) {
    term *= x/(a+n);
    sum  += term;
    if(fabs(term/sum) < GSL_DBL_EPSILON) break;
  }

  result->val  = D.val * sum;
  result->err  = D.err * fabs(sum);
  result->err += (1.0 + n) * GSL_DBL_EPSILON * fabs(result->val);

  if(n == nmax)
    GSL_ERROR ("error", GSL_EMAXITER);
  else
    return stat_D;
}

/* Q large x asymptotic
 */
static
int
gamma_inc_Q_large_x(const double a, const double x, gsl_sf_result * result)
{
  const int nmax = 5000;

  gsl_sf_result D;
  const int stat_D = gamma_inc_D(a, x, &D);

  double sum  = 1.0;
  double term = 1.0;
  double last = 1.0;
  int n;
  for(n=1; n<nmax; n++) {
    term *= (a-n)/x;
    if(fabs(term/last) > 1.0) break;
    if(fabs(term/sum)  < GSL_DBL_EPSILON) break;
    sum  += term;
    last  = term;
  }

  result->val  = D.val * (a/x) * sum;
  result->err  = D.err * fabs((a/x) * sum);
  result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);

  if(n == nmax)
    GSL_ERROR ("error in large x asymptotic", GSL_EMAXITER);
  else
    return stat_D;
}

/* Uniform asymptotic for x near a, a and x large.
 * See [Temme, p. 285]
 */
static
int
gamma_inc_Q_asymp_unif(const double a, const double x, gsl_sf_result * result)
{
  const double rta = sqrt(a);
  const double eps = (x-a)/a;

  gsl_sf_result ln_term;
  const int stat_ln = gsl_sf_log_1plusx_mx_e(eps, &ln_term);  /* log(1+eps) - eps */
  const double eta  = GSL_SIGN(eps) * sqrt(-2.0*ln_term.val);

  gsl_sf_result erfc;

  double R;
  double c0, c1;

  /* This used to say erfc(eta*M_SQRT2*rta), which is wrong.
   * The sqrt(2) is in the denominator. Oops.
   * Fixed: [GJ] Mon Nov 15 13:25:32 MST 2004
   */
  gsl_sf_erfc_e(eta*rta/M_SQRT2, &erfc);

  if(fabs(eps) < GSL_ROOT5_DBL_EPSILON) {
    c0 = -1.0/3.0 + eps*(1.0/12.0 - eps*(23.0/540.0 - eps*(353.0/12960.0 - eps*589.0/30240.0)));
    c1 = -1.0/540.0 - eps/288.0;
  }
  else {
    const double rt_term = sqrt(-2.0 * ln_term.val/(eps*eps));
    const double lam = x/a;
    c0 = (1.0 - 1.0/rt_term)/eps;
    c1 = -(eta*eta*eta * (lam*lam + 10.0*lam + 1.0) - 12.0 * eps*eps*eps) / (12.0 * eta*eta*eta*eps*eps*eps);
  }

  R = exp(-0.5*a*eta*eta)/(M_SQRT2*M_SQRTPI*rta) * (c0 + c1/a);

  result->val  = 0.5 * erfc.val + R;
  result->err  = GSL_DBL_EPSILON * fabs(R * 0.5 * a*eta*eta) + 0.5 * erfc.err;
  result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);

  return stat_ln;
}

/* Continued fraction which occurs in evaluation
 * of Q(a,x) or Gamma(a,x).
 *
 *              1   (1-a)/x  1/x  (2-a)/x   2/x  (3-a)/x
 *   F(a,x) =  ---- ------- ----- -------- ----- -------- ...
 *             1 +   1 +     1 +   1 +      1 +   1 +
 *
 * Hans E. Plesser, 2002-01-22 (hans dot plesser at itf dot nlh dot no).
 *
 * Split out from gamma_inc_Q_CF() by GJ [Tue Apr  1 13:16:41 MST 2003].
 * See gamma_inc_Q_CF() below.
 *
 */
static int
gamma_inc_F_CF(const double a, const double x, gsl_sf_result * result)
{
  const int    nmax  =  5000;
  const double small =  gsl_pow_3 (GSL_DBL_EPSILON);

  double hn = 1.0;           /* convergent */
  double Cn = 1.0 / small;
  double Dn = 1.0;
  int n;

  /* n == 1 has a_1, b_1, b_0 independent of a,x,
     so that has been done by hand                */
  for ( n = 2 ; n < nmax ; n++ )
  {
    double an;
    double delta;

    if(GSL_IS_ODD(n))
      an = 0.5*(n-1)/x;
    else
      an = (0.5*n-a)/x;

    Dn = 1.0 + an * Dn;
    if ( fabs(Dn) < small )
      Dn = small;
    Cn = 1.0 + an/Cn;
    if ( fabs(Cn) < small )
      Cn = small;
    Dn = 1.0 / Dn;
    delta = Cn * Dn;
    hn *= delta;
    if(fabs(delta-1.0) < GSL_DBL_EPSILON) break;
  }

  result->val = hn;
  result->err = 2.0*GSL_DBL_EPSILON * fabs(hn);
  result->err += GSL_DBL_EPSILON * (2.0 + 0.5*n) * fabs(result->val);

  if(n == nmax)
    GSL_ERROR ("error in CF for F(a,x)", GSL_EMAXITER);
  else
    return GSL_SUCCESS;
}

/* Continued fraction for Q.
 *
 * Q(a,x) = D(a,x) a/x F(a,x)
 *
 * Hans E. Plesser, 2002-01-22 (hans dot plesser at itf dot nlh dot no):
 *
 * Since the Gautschi equivalent series method for CF evaluation may lead
 * to singularities, I have replaced it with the modified Lentz algorithm
 * given in
 *
 * I J Thompson and A R Barnett
 * Coulomb and Bessel Functions of Complex Arguments and Order
 * J Computational Physics 64:490-509 (1986)
 *
 * In consequence, gamma_inc_Q_CF_protected() is now obsolete and has been
 * removed.
 *
 * Identification of terms between the above equation for F(a, x) and
 * the first equation in the appendix of Thompson&Barnett is as follows:
 *
 *    b_0 = 0, b_n = 1 for all n > 0
 *
 *    a_1 = 1
 *    a_n = (n/2-a)/x    for n even
 *    a_n = (n-1)/(2x)   for n odd
 *
 */
static
int
gamma_inc_Q_CF(const double a, const double x, gsl_sf_result * result)
{
  gsl_sf_result D;
  gsl_sf_result F;
  const int stat_D = gamma_inc_D(a, x, &D);
  const int stat_F = gamma_inc_F_CF(a, x, &F);

  result->val  = D.val * (a/x) * F.val;
  result->err  = D.err * fabs((a/x) * F.val) + fabs(D.val * a/x * F.err);

  return GSL_ERROR_SELECT_2(stat_F, stat_D);
}

/* Useful for small a and x. Handles the subtraction analytically.
 */
static
int
gamma_inc_Q_series(const double a, const double x, gsl_sf_result * result)
{
  double term1;  /* 1 - x^a/Gamma(a+1) */
  double sum;    /* 1 + (a+1)/(a+2)(-x)/2! + (a+1)/(a+3)(-x)^2/3! + ... */
  int stat_sum;
  double term2;  /* a temporary variable used at the end */

  {
    /* Evaluate series for 1 - x^a/Gamma(a+1), small a
     */
    const double pg21 = -2.404113806319188570799476;  /* PolyGamma[2,1] */
    const double lnx  = log(x);
    const double el   = M_EULER+lnx;
    const double c1 = -el;
    const double c2 = M_PI*M_PI/12.0 - 0.5*el*el;
    const double c3 = el*(M_PI*M_PI/12.0 - el*el/6.0) + pg21/6.0;
    const double c4 = -0.04166666666666666667
                       * (-1.758243446661483480 + lnx)
                       * (-0.764428657272716373 + lnx)
                       * ( 0.723980571623507657 + lnx)
                       * ( 4.107554191916823640 + lnx);
    const double c5 = -0.0083333333333333333
                       * (-2.06563396085715900 + lnx)
                       * (-1.28459889470864700 + lnx)
                       * (-0.27583535756454143 + lnx)
                       * ( 1.33677371336239618 + lnx)
                       * ( 5.17537282427561550 + lnx);
    const double c6 = -0.0013888888888888889
                       * (-2.30814336454783200 + lnx)
                       * (-1.65846557706987300 + lnx)
                       * (-0.88768082560020400 + lnx)
                       * ( 0.17043847751371778 + lnx)
                       * ( 1.92135970115863890 + lnx)
                       * ( 6.22578557795474900 + lnx);
    const double c7 = -0.00019841269841269841
                       * (-2.5078657901291800 + lnx)
                       * (-1.9478900888958200 + lnx)
                       * (-1.3194837322612730 + lnx)
                       * (-0.5281322700249279 + lnx)
                       * ( 0.5913834939078759 + lnx)
                       * ( 2.4876819633378140 + lnx)
                       * ( 7.2648160783762400 + lnx);
    const double c8 = -0.00002480158730158730
                       * (-2.677341544966400 + lnx)
                       * (-2.182810448271700 + lnx)
                       * (-1.649350342277400 + lnx)
                       * (-1.014099048290790 + lnx)
                       * (-0.191366955370652 + lnx)
                       * ( 0.995403817918724 + lnx)
                       * ( 3.041323283529310 + lnx)
                       * ( 8.295966556941250 + lnx);
    const double c9 = -2.75573192239859e-6
                       * (-2.8243487670469080 + lnx)
                       * (-2.3798494322701120 + lnx)
                       * (-1.9143674728689960 + lnx)
                       * (-1.3814529102920370 + lnx)
                       * (-0.7294312810261694 + lnx)
                       * ( 0.1299079285269565 + lnx)
                       * ( 1.3873333251885240 + lnx)
                       * ( 3.5857258865210760 + lnx)
                       * ( 9.3214237073814600 + lnx);
    const double c10 = -2.75573192239859e-7
                       * (-2.9540329644556910 + lnx)
                       * (-2.5491366926991850 + lnx)
                       * (-2.1348279229279880 + lnx)
                       * (-1.6741881076349450 + lnx)
                       * (-1.1325949616098420 + lnx)
                       * (-0.4590034650618494 + lnx)
                       * ( 0.4399352987435699 + lnx)
                       * ( 1.7702236517651670 + lnx)
                       * ( 4.1231539047474080 + lnx)
                       * ( 10.342627908148680 + lnx);

    term1 = a*(c1+a*(c2+a*(c3+a*(c4+a*(c5+a*(c6+a*(c7+a*(c8+a*(c9+a*c10)))))))));
  }

  {
    /* Evaluate the sum.
     */
    const int nmax = 5000;
    double t = 1.0;
    int n;
    sum = 1.0;

    for(n=1; n<nmax; n++) {
      t *= -x/(n+1.0);
      sum += (a+1.0)/(a+n+1.0)*t;
      if(fabs(t/sum) < GSL_DBL_EPSILON) break;
    }

    if(n == nmax)
      stat_sum = GSL_EMAXITER;
    else
      stat_sum = GSL_SUCCESS;
  }

  term2 = (1.0 - term1) * a/(a+1.0) * x * sum;
  result->val  = term1 + term2;
  result->err  = GSL_DBL_EPSILON * (fabs(term1) + 2.0*fabs(term2));
  result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  return stat_sum;
}

/* series for small a and x, but not defined for a == 0 */
static int
gamma_inc_series(double a, double x, gsl_sf_result * result)
{
  gsl_sf_result Q;
  gsl_sf_result G;
  const int stat_Q = gamma_inc_Q_series(a, x, &Q);
  const int stat_G = gsl_sf_gamma_e(a, &G);
  result->val = Q.val * G.val;
  result->err = fabs(Q.val * G.err) + fabs(Q.err * G.val);
  result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);

  return GSL_ERROR_SELECT_2(stat_Q, stat_G);
}

static int
gamma_inc_a_gt_0(double a, double x, gsl_sf_result * result)
{
  /* x > 0 and a > 0; use result for Q */
  gsl_sf_result Q;
  gsl_sf_result G;
  const int stat_Q = gsl_sf_gamma_inc_Q_e(a, x, &Q);
  const int stat_G = gsl_sf_gamma_e(a, &G);

  result->val = G.val * Q.val;
  result->err = fabs(G.val * Q.err) + fabs(G.err * Q.val);
  result->err += 2.0*GSL_DBL_EPSILON * fabs(result->val);

  return GSL_ERROR_SELECT_2(stat_G, stat_Q);
}

static int
gamma_inc_CF(double a, double x, gsl_sf_result * result)
{
  gsl_sf_result F;
  gsl_sf_result pre;
  const int stat_F = gamma_inc_F_CF(a, x, &F);
  const int stat_E = gsl_sf_exp_e((a-1.0)*log(x) - x, &pre);

  result->val = F.val * pre.val;
  result->err = fabs(F.err * pre.val) + fabs(F.val * pre.err);
  result->err += (2.0 + fabs(a)) * GSL_DBL_EPSILON * fabs(result->val);

  return GSL_ERROR_SELECT_2(stat_F, stat_E);
}

/* evaluate Gamma(0,x), x > 0 */
#define GAMMA_INC_A_0(x, result) gsl_sf_expint_E1_e(x, result)

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

int
gsl_sf_gamma_inc_Q_e(const double a, const double x, gsl_sf_result * result)
{
  if(a < 0.0 || x < 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x == 0.0) {
    result->val = 1.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(a == 0.0)
  {
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(x <= 0.5*a) {
    /* If the series is quick, do that. It is
     * robust and simple.
     */
    gsl_sf_result P;
    int stat_P = gamma_inc_P_series(a, x, &P);
    result->val  = 1.0 - P.val;
    result->err  = P.err;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return stat_P;
  }
  else if(a >= 1.0e+06 && (x-a)*(x-a) < a) {
    /* Then try the difficult asymptotic regime.
     * This is the only way to do this region.
     */
    return gamma_inc_Q_asymp_unif(a, x, result);
  }
  else if(a < 0.2 && x < 5.0) {
    /* Cancellations at small a must be handled
     * analytically; x should not be too big
     * either since the series terms grow
     * with x and log(x).
     */
    return gamma_inc_Q_series(a, x, result);
  }
  else if(a <= x) {
    if(x <= 1.0e+06) {
      /* Continued fraction is excellent for x >~ a.
       * We do not let x be too large when x > a since
       * it is somewhat pointless to try this there;
       * the function is rapidly decreasing for
       * x large and x > a, and it will just
       * underflow in that region anyway. We
       * catch that case in the standard
       * large-x method.
       */
      return gamma_inc_Q_CF(a, x, result);
    }
    else {
      return gamma_inc_Q_large_x(a, x, result);
    }
  }
  else {
    if(x > a - sqrt(a)) {
      /* Continued fraction again. The convergence
       * is a little slower here, but that is fine.
       * We have to trade that off against the slow
       * convergence of the series, which is the
       * only other option.
       */
      return gamma_inc_Q_CF(a, x, result);
    }
    else {
      gsl_sf_result P;
      int stat_P = gamma_inc_P_series(a, x, &P);
      result->val  = 1.0 - P.val;
      result->err  = P.err;
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
      return stat_P;
    }
  }
}

int
gsl_sf_gamma_inc_P_e(const double a, const double x, gsl_sf_result * result)
{
  if(a <= 0.0 || x < 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x == 0.0) {
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(x < 20.0 || x < 0.5*a) {
    /* Do the easy series cases. Robust and quick.
     */
    return gamma_inc_P_series(a, x, result);
  }
  else if(a > 1.0e+06 && (x-a)*(x-a) < a) {
    /* Crossover region. Note that Q and P are
     * roughly the same order of magnitude here,
     * so the subtraction is stable.
     */
    gsl_sf_result Q;
    int stat_Q = gamma_inc_Q_asymp_unif(a, x, &Q);
    result->val  = 1.0 - Q.val;
    result->err  = Q.err;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return stat_Q;
  }
  else if(a <= x) {
    /* Q <~ P in this area, so the
     * subtractions are stable.
     */
    gsl_sf_result Q;
    int stat_Q;
    if(a > 0.2*x) {
      stat_Q = gamma_inc_Q_CF(a, x, &Q);
    }
    else {
      stat_Q = gamma_inc_Q_large_x(a, x, &Q);
    }
    result->val  = 1.0 - Q.val;
    result->err  = Q.err;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return stat_Q;
  }
  else {
    if((x-a)*(x-a) < a) {
      /* This condition is meant to insure
       * that Q is not very close to 1,
       * so the subtraction is stable.
       */
      gsl_sf_result Q;
      int stat_Q = gamma_inc_Q_CF(a, x, &Q);
      result->val  = 1.0 - Q.val;
      result->err  = Q.err;
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
      return stat_Q;
    }
    else {
      return gamma_inc_P_series(a, x, result);
    }
  }
}

int
gsl_sf_gamma_inc_e(const double a, const double x, gsl_sf_result * result)
{
  if(x < 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x == 0.0) {
    return gsl_sf_gamma_e(a, result);
  }
  else if(a == 0.0)
  {
    return GAMMA_INC_A_0(x, result);
  }
  else if(a > 0.0)
  {
    return gamma_inc_a_gt_0(a, x, result);
  }
  else if(x > 0.25)
  {
    /* continued fraction seems to fail for x too small; otherwise
       it is ok, independent of the value of |x/a|, because of the
       non-oscillation in the expansion, i.e. the CF is
       un-conditionally convergent for a < 0 and x > 0
     */
    return gamma_inc_CF(a, x, result);
  }
  else if(fabs(a) < 0.5)
  {
    return gamma_inc_series(a, x, result);
  }
  else
  {
    /* a = fa + da; da >= 0 */
    const double fa = floor(a);
    const double da = a - fa;

    gsl_sf_result g_da;
    const int stat_g_da = ( da > 0.0 ? gamma_inc_a_gt_0(da, x, &g_da)
                                     : GAMMA_INC_A_0(x, &g_da));

    double alpha = da;
    double gax = g_da.val;

    /* Gamma(alpha-1,x) = 1/(alpha-1) (Gamma(a,x) - x^(alpha-1) e^-x) */
    do
    {
      const double shift = exp(-x + (alpha-1.0)*log(x));
      gax = (gax - shift) / (alpha - 1.0);
      alpha -= 1.0;
    } while(alpha > a);

    result->val = gax;
    result->err = 2.0*(1.0 + fabs(a))*GSL_DBL_EPSILON*fabs(gax);
    return stat_g_da;
  }

}

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

#include "eval.h"

double gsl_sf_gamma_inc_P(const double a, const double x)
{
  EVAL_RESULT(gsl_sf_gamma_inc_P_e(a, x, &result));
}

double gsl_sf_gamma_inc_Q(const double a, const double x)
{
  EVAL_RESULT(gsl_sf_gamma_inc_Q_e(a, x, &result));
}

double gsl_sf_gamma_inc(const double a, const double x)
{
   EVAL_RESULT(gsl_sf_gamma_inc_e(a, x, &result));
}


text:
/*
ODE: a program to get optime Runge-Kutta and multi-steps methods.

Copyright 2011-2019, Javier Burguete Tolosa.

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 Javier Burguete Tolosa ``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 Javier Burguete Tolosa 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 optimize.c
 * \brief Source file with optimize functions.
 * \author Javier Burguete Tolosa.
 * \copyright Copyright 2011-2019.
 */
#define _GNU_SOURCE
#include <stdio.h>
#include <string.h>
#include <float.h>
#include <math.h>
#include <libxml/parser.h>
#include <glib.h>
#include <libintl.h>
#include <gsl/gsl_rng.h>
#if HAVE_MPI
#include <mpi.h>
#endif
#include "config.h"
#include "utils.h"
#include "optimize.h"

#define DEBUG_OPTIMIZE 0        ///< macro to debug.

GMutex mutex[1];                ///< GMutex struct.
FILE *file_variables = NULL;    ///< random variables file.
int rank;                       ///< MPI rank.
int nnodes;                     ///< MPI nodes number.
unsigned nthreads;              ///< threads number.

/**
 * Function to print the random variables on a file.
 */
void
optimize_print_random (Optimize * optimize,     ///< Optimize struct.
                       FILE * file)     ///< file.
{
  unsigned int i, n;
  n = optimize->nfree;
  for (i = 0; i < n; ++i)
    fprintf (file, "o%d:%.19Le;\n", i, optimize->value_optimal[i]);
  for (i = 0; i < n; ++i)
    fprintf (file, "m%d:%.19Le;\n", i, optimize->minimum[i]);
  for (i = 0; i < n; ++i)
    fprintf (file, "i%d:%.19Le;\n", i, optimize->interval[i]);
}

/**
 * Function to perform every optimization step.
 */
void
optimize_step (Optimize * optimize)     ///< Optimize struct.
{
  long double *is, *vo, *vo2, *random;
  long double o, o2, v, f;
  unsigned long long int ii, nrandom;
  unsigned int i, j, k, n, nfree;

#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_step: start\n");
#endif

  // save optimal values
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_step: save optimal values\n");
#endif
  nfree = optimize->nfree;
  o2 = *optimize->optimal;
  vo = (long double *) alloca (nfree * sizeof (long double));
  vo2 = (long double *) alloca (nfree * sizeof (long double));
  memcpy (vo, optimize->value_optimal, nfree * sizeof (long double));

  // optimization algorithm sampling
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_step: optimization algorithm sampling\n");
  fprintf (stderr, "optimize_step: nsimulations=%Lu\n", optimize->nsimulations);
#endif
  random = optimize->random_data;
  ii = optimize->nsimulations * (rank * nthreads + optimize->thread)
    / (nnodes * nthreads);
  nrandom = optimize->nsimulations * (rank * nthreads + optimize->thread + 1)
    / (nnodes * nthreads);
  for (; ii < nrandom; ++ii)
    {

      // random freedom degrees
#if DEBUG_OPTIMIZE
      fprintf (stderr, "optimize_step: random freedom degrees\n");
      fprintf (stderr, "optimize_step: simulation=%Lu\n", ii);
#endif
      optimize_generate_freedom (optimize, ii);

      // method coefficients
#if DEBUG_OPTIMIZE
      fprintf (stderr, "optimize_step: method coefficients\n");
#endif
      if (!optimize->method (optimize))
        o = INFINITY;
      else
        o = optimize->objective (optimize);
      if (o < o2)
        {
          o2 = o;
          memcpy (vo, random, nfree * sizeof (long double));
        }
      if (file_variables)
        {
          g_mutex_lock (mutex);
          print_variables (random, nfree, file_variables);
          fprintf (file_variables, "%.19Le\n", o);
          g_mutex_unlock (mutex);
        }
    }

  // array of intervals to climb around the optimal
#if DEBUG_OPTIMIZE
  fprintf (stderr,
           "optimize_step: array of intervals to climb around the optimal\n");
#endif
  is = (long double *) alloca (nfree * sizeof (long double));
  for (j = 0; j < nfree; ++j)
    is[j] = optimize->interval0[j] * optimize->climbing_factor;

  // hill climbing algorithm bucle
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_step: hill climbing algorithm bucle\n");
#endif
  memcpy (vo2, vo, nfree * sizeof (long double));
  n = optimize->nclimbings;
  for (i = 0; i < n; ++i)
    {
      memcpy (random, vo, nfree * sizeof (long double));
      for (j = k = 0; j < nfree; ++j)
        {
          v = vo[j];
          random[j] = v + is[j];
          if (!optimize->method (optimize))
            o = INFINITY;
          else
            o = optimize->objective (optimize);
          if (o < o2)
            {
              k = 1;
              o2 = o;
              memcpy (vo2, random, nfree * sizeof (long double));
            }
          if (file_variables)
            {
              g_mutex_lock (mutex);
              print_variables (random, nfree, file_variables);
              fprintf (file_variables, "%.19Le\n", o);
              g_mutex_unlock (mutex);
            }
          random[j] = fmaxl (0.L, v - is[j]);
          if (!optimize->method (optimize))
            o = INFINITY;
          else
            o = optimize->objective (optimize);
          if (o < o2)
            {
              k = 1;
              o2 = o;
              memcpy (vo2, random, nfree * sizeof (long double));
            }
          if (file_variables)
            {
              g_mutex_lock (mutex);
              print_variables (random, nfree, file_variables);
              fprintf (file_variables, "%.19Le\n", o);
              g_mutex_unlock (mutex);
            }
          random[j] = v;
        }

      // update optimal values and increase or reduce intervals if converging or
      // not
      if (!k)
        f = 0.5L;
      else
        {
          f = 1.2L;
          memcpy (vo, vo2, nfree * sizeof (long double));
        }
      for (j = 0; j < nfree; ++j)
        is[j] *= f;
    }

  // update optimal values
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_step: update optimal values\n");
#endif
  if (o2 < *optimize->optimal)
    {
      g_mutex_lock (mutex);
      *optimize->optimal = o2;
      memcpy (optimize->value_optimal, vo2, nfree * sizeof (long double));
      g_mutex_unlock (mutex);
    }

#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_step: end\n");
#endif
}

/**
 * Function to init required variables on an Optimize struct data.
 */
void
optimize_init (Optimize * optimize,     ///< Optimize struct.
               gsl_rng * rng,   ///< GSL pseudo-random number generator struct.
               unsigned int thread)     ///< thread number.
{
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_init: start\n");
#endif
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_init: nsimulations=%Lu nfree=%u size=%u\n",
           optimize->nsimulations, optimize->nfree, optimize->size);
#endif
  optimize->random_data
    = (long double *) g_slice_alloc (optimize->nfree * sizeof (long double));
  optimize->coefficient
    = (long double *) g_slice_alloc (optimize->size * sizeof (long double));
  optimize->minimum
    = (long double *) g_slice_alloc (optimize->nfree * sizeof (long double));
  optimize->interval
    = (long double *) g_slice_alloc (optimize->nfree * sizeof (long double));
  memcpy (optimize->minimum, optimize->minimum0,
          optimize->nfree * sizeof (long double));
  memcpy (optimize->interval, optimize->interval0,
          optimize->nfree * sizeof (long double));
  optimize->rng = rng;
  optimize->thread = thread;
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_init: end\n");
#endif
}

/**
 * Function to free the memory allocated by an Optimize struct.
 */
void
optimize_delete (Optimize * optimize)   ///< Optimize struct.
{
  g_slice_free1 (optimize->nfree * sizeof (long double), optimize->interval);
  g_slice_free1 (optimize->nfree * sizeof (long double), optimize->minimum);
  g_slice_free1 (optimize->size * sizeof (long double), optimize->coefficient);
  g_slice_free1 (optimize->nfree * sizeof (long double), optimize->random_data);
}

/**
 * Function to do the optimization bucle.
 */
void
optimize_bucle (Optimize * optimize)    ///< Optimize struct.
{
  GThread *thread[nthreads];
#if HAVE_MPI
  long double *vo;
  MPI_Status status;
#endif
  unsigned int i, j, nfree;

#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_bucle: start\n");
  fprintf (stderr, "optimize_bucle: nfree=%u\n", optimize->nfree);
#endif

  // Allocate local array of optimal values
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_bucle: allocate local array of optimal values\n");
#endif
  nfree = optimize->nfree;
#if HAVE_MPI
  vo = (long double *) alloca ((1 + nfree) * sizeof (long double));
  printf ("Rank=%d NNodes=%d\n", rank, nnodes);
#endif

  // Init some parameters
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_bucle: init some parameters\n");
#endif
  *optimize->optimal = INFINITY;
  for (i = 0; i < nfree; ++i)
    optimize->value_optimal[i]
      = optimize->minimum[i] + 0.5L * optimize->interval[i];

  // Iterate
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_bucle: iterate\n");
#endif
  for (i = 0; i < optimize->niterations; ++i)
    {

      // Optimization step parallelized for every node by GThreads
      if (nthreads > 1)
        {
          for (j = 0; j < nthreads; ++j)
            thread[j]
              = g_thread_new (NULL,
                              (GThreadFunc) (void (*)(void)) optimize_step,
                              (void *) (optimize + j));
          for (j = 0; j < nthreads; ++j)
            g_thread_join (thread[j]);
        }
      else
        optimize_step (optimize);

#if HAVE_MPI
      if (rank > 0)
        {

          // Secondary nodes send the optimal coefficients to the master node
          vo[0] = *optimize->optimal;
          memcpy (vo + 1, optimize->value_optimal,
                  nfree * sizeof (long double));
          MPI_Send (vo, 1 + nfree, MPI_LONG_DOUBLE, 0, 1, MPI_COMM_WORLD);

          // Secondary nodes receive the optimal coefficients
          MPI_Recv (vo, 1 + nfree, MPI_LONG_DOUBLE, 0, 1, MPI_COMM_WORLD,
                    &status);
          *optimize->optimal = vo[0];
          memcpy (optimize->value_optimal, vo + 1,
                  nfree * sizeof (long double));
        }
      else
        {
          printf ("rank=%d optimal=%.19Le\n", rank, *optimize->optimal);

          for (j = 1; j < nnodes; ++j)
            {

              // Master node receives the optimal coefficients obtained by
              // secondary nodes
              MPI_Recv (vo, 1 + nfree, MPI_LONG_DOUBLE, j, 1, MPI_COMM_WORLD,
                        &status);

              // Master node selects the optimal coefficients
              if (vo[0] < *optimize->optimal)
                {
                  *optimize->optimal = vo[0];
                  memcpy (optimize->value_optimal, vo + 1,
                          nfree * sizeof (long double));
                }
            }

          // Master node sends the optimal coefficients to secondary nodes
          vo[0] = *optimize->optimal;
          memcpy (vo + 1, optimize->value_optimal,
                  nfree * sizeof (long double));
          for (j = 1; j < nnodes; ++j)
            MPI_Send (vo, 1 + nfree, MPI_LONG_DOUBLE, j, 1, MPI_COMM_WORLD);
        }

#endif

      // Print the optimal coefficients
#if DEBUG_OPTIMIZE
      optimize_print_random (optimize, stderr);
      fprintf (stderr, "optimal=%.19Le\n", *optimize->optimal);
#endif

      // Updating coefficient intervals to converge
      optimize_converge (optimize);

      // Iterate
#if HAVE_MPI
      printf ("Rank %u\n", rank);
#endif
      printf ("Iteration %u Optimal %.19Le\n", i + 1, *optimize->optimal);
    }

#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_bucle: end\n");
#endif
}

/**
 * Function to create an Optimize struct data.
 */
void
optimize_create (Optimize * optimize,   ///< Optimize struct.
                 long double *optimal,
                 ///< pointer to the optimal objective function value.
                 long double *value_optimal)
                 ///< array of optimal freedom degree values.
{
  unsigned long long int nsimulations;
  unsigned int i, nfree;
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_create: start\n");
#endif
  optimize->optimal = optimal;
  optimize->value_optimal = value_optimal;
  nfree = optimize->nfree;
  optimize->nsimulations = nsimulations = optimize->nvariable;
  for (i = 1; i < nfree; ++i)
    optimize->nsimulations *= nsimulations;
  optimize->nclimbings *= nfree;
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_create nsimulations=%Lu nclimbings=%u nfree=%u\n",
           optimize->nsimulations, optimize->nclimbings, nfree);
  fprintf (stderr, "optimize_create: end\n");
#endif
}

/**
 * Function to read the Optimize struct data on a XML node.
 *
 * \return 1 on success, 0 on error.
 */
int
optimize_read (Optimize * optimize,     ///< Optimize struct.
               xmlNode * node)  ///< XML node.
{
  int code;
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_read: start\n");
#endif
  optimize->nvariable = xml_node_get_uint (node, XML_NSIMULATIONS, &code);
  if (code || !optimize->nvariable)
    {
      error_message = g_strdup (_("Bad simulations number"));
      goto exit_on_error;
    }
  optimize->nclimbings
    = xml_node_get_uint_with_default (node, XML_NCLIMBINGS, 0, &code);
  if (code)
    {
      error_message = g_strdup (_("Bad hill climbings number"));
      goto exit_on_error;
    }
  optimize->niterations = xml_node_get_uint (node, XML_NITERATIONS, &code);
  if (code || !optimize->niterations)
    {
      error_message = g_strdup (_("Bad iterations number"));
      goto exit_on_error;
    }
  optimize->convergence_factor
    = xml_node_get_float (node, XML_CONVERGENCE_FACTOR, &code);
  if (code || optimize->convergence_factor < LDBL_EPSILON)
    {
      error_message = g_strdup (_("Bad convergence factor"));
      goto exit_on_error;
    }
  optimize->climbing_factor
    = xml_node_get_float (node, XML_CLIMBING_FACTOR, &code);
  if (code || optimize->climbing_factor < LDBL_EPSILON)
    {
      error_message = g_strdup (_("Bad climging factor"));
      goto exit_on_error;
    }
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_read: end\n");
#endif
  return 1;

exit_on_error:
#if DEBUG_OPTIMIZE
  fprintf (stderr, "optimize_read: end\n");
#endif
  return 0;
}

max_stars_repo_name: jburguete/ode

max_issues_repo_name: jburguete/ode

max_forks_repo_name: jburguete/ode


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

#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_movstat.h>

/* compute moving QQR by explicitely constructing window and computing QQR */
int
slow_movqqr(const gsl_movstat_end_t etype, const double q, const gsl_vector * x,
            gsl_vector * y, const int H, const int J)
{
  const size_t n = x->size;
  const int K = H + J + 1;
  double *window = malloc(K * sizeof(double));
  size_t i;

  for (i = 0; i < n; ++i)
    {
      size_t wsize = gsl_movstat_fill(etype, x, i, H, J, window);
      double quant1, quant2;

      gsl_sort(window, 1, wsize);
      quant1 = gsl_stats_quantile_from_sorted_data(window, 1, wsize, q);
      quant2 = gsl_stats_quantile_from_sorted_data(window, 1, wsize, 1.0 - q);

      gsl_vector_set(y, i, quant2 - quant1);
    }

  free(window);

  return GSL_SUCCESS;
}

static double
func_qqr(const size_t n, double x[], void * params)
{
  double q = *(double *) params;
  double quant1, quant2;

  (void) params;

  gsl_sort(x, 1, n);
  quant1 = gsl_stats_quantile_from_sorted_data(x, 1, n, q);
  quant2 = gsl_stats_quantile_from_sorted_data(x, 1, n, 1.0 - q);

  return (quant2 - quant1);
}

static void
test_qqr_proc(const double tol, const double q, const size_t n, const size_t H, const size_t J,
              const gsl_movstat_end_t etype, gsl_rng * rng_p)
{
  gsl_movstat_workspace * w = gsl_movstat_alloc2(H, J);
  gsl_vector * x = gsl_vector_alloc(n);
  gsl_vector * y = gsl_vector_alloc(n);
  gsl_vector * z = gsl_vector_alloc(n);
  gsl_movstat_function F;
  char buf[2048];

  F.function = func_qqr;
  F.params = (void *) &q;

  random_vector(x, rng_p);

  /* y = QQR(x) with slow brute force method */
  slow_movqqr(etype, q, x, y, H, J);

  /* y = QQR(x) with fast method */
  gsl_movstat_qqr(etype, x, q, z, w);

  /* test y = z */
  sprintf(buf, "n=%zu H=%zu J=%zu endtype=%u QQR random", n, H, J, etype);
  compare_vectors(tol, z, y, buf);

  /* z = QQR(x) in-place */
  gsl_vector_memcpy(z, x);
  gsl_movstat_qqr(etype, z, q, z, w);

  sprintf(buf, "n=%zu H=%zu J=%zu endtype=%u QQR random in-place", n, H, J, etype);
  compare_vectors(tol, z, y, buf);

  /* z = QQR(x) with user-defined function */
  gsl_movstat_apply(etype, &F, x, z, w);

  sprintf(buf, "n=%zu H=%zu J=%zu endtype=%u QQR user", n, H, J, etype);
  compare_vectors(tol, z, y, buf);

  gsl_movstat_free(w);
  gsl_vector_free(x);
  gsl_vector_free(y);
  gsl_vector_free(z);
}

static void
test_qqr(gsl_rng * rng_p)
{
  const double eps = 1.0e-10;

  test_qqr_proc(eps, 0.1, 100, 0, 0, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_qqr_proc(eps, 0.1, 1000, 3, 3, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_qqr_proc(eps, 0.2, 1000, 0, 5, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_qqr_proc(eps, 0.3, 1000, 5, 0, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_qqr_proc(eps, 0.25, 2000, 10, 5, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_qqr_proc(eps, 0.4, 2000, 5, 10, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_qqr_proc(eps, 0.25, 20, 50, 50, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_qqr_proc(eps, 0.25, 20, 10, 50, GSL_MOVSTAT_END_PADZERO, rng_p);
  test_qqr_proc(eps, 0.25, 20, 50, 10, GSL_MOVSTAT_END_PADZERO, rng_p);

  test_qqr_proc(eps, 0.1, 100, 0, 0, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_qqr_proc(eps, 0.1, 1000, 3, 3, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_qqr_proc(eps, 0.2, 1000, 0, 5, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_qqr_proc(eps, 0.3, 1000, 5, 0, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_qqr_proc(eps, 0.25, 2000, 10, 5, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_qqr_proc(eps, 0.4, 2000, 5, 10, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_qqr_proc(eps, 0.25, 20, 50, 50, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_qqr_proc(eps, 0.25, 20, 10, 50, GSL_MOVSTAT_END_PADVALUE, rng_p);
  test_qqr_proc(eps, 0.25, 20, 50, 10, GSL_MOVSTAT_END_PADVALUE, rng_p);

  test_qqr_proc(eps, 0.1, 100, 0, 0, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_qqr_proc(eps, 0.1, 1000, 3, 3, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_qqr_proc(eps, 0.2, 1000, 0, 5, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_qqr_proc(eps, 0.3, 1000, 5, 0, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_qqr_proc(eps, 0.25, 2000, 10, 5, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_qqr_proc(eps, 0.4, 2000, 5, 10, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_qqr_proc(eps, 0.25, 20, 50, 50, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_qqr_proc(eps, 0.25, 20, 10, 50, GSL_MOVSTAT_END_TRUNCATE, rng_p);
  test_qqr_proc(eps, 0.25, 20, 50, 10, GSL_MOVSTAT_END_TRUNCATE, rng_p);
}