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/GraphicsAttributes.h>
#include <terminal/Hyperlink.h>
#include <terminal/primitives.h>

#include <crispy/BufferObject.h>
#include <crispy/Comparison.h>
#include <crispy/assert.h>

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

#include <iterator>
#include <sstream>
#include <string>
#include <vector>

namespace terminal
{

enum class LineFlags : uint8_t
{
    None = 0x0000,
    Wrappable = 0x0001,
    Wrapped = 0x0002,
    Marked = 0x0004,
    // TODO: DoubleWidth  = 0x0010,
    // TODO: DoubleHeight = 0x0020,
};

// clang-format off
template <typename, bool> struct OptionalProperty;
template <typename T> struct OptionalProperty<T, false> {};
template <typename T> struct OptionalProperty<T, true> { T value; };
// clang-format on

/**
 * Line storage with call columns sharing the same SGR attributes.
 */
struct TrivialLineBuffer
{
    ColumnCount displayWidth;
    GraphicsAttributes textAttributes;
    GraphicsAttributes fillAttributes = textAttributes;
    HyperlinkId hyperlink {};

    ColumnCount usedColumns {};
    crispy::BufferFragment text {};

    void reset(GraphicsAttributes _attributes) noexcept
    {
        textAttributes = _attributes;
        fillAttributes = _attributes;
        hyperlink = {};
        usedColumns = {};
        text.reset();
    }
};

template <typename Cell>
using InflatedLineBuffer = std::vector<Cell>;

/// Unpacks a TrivialLineBuffer into an InflatedLineBuffer<Cell>.
template <typename Cell>
InflatedLineBuffer<Cell> inflate(TrivialLineBuffer const& input);

template <typename Cell>
using LineStorage = std::variant<TrivialLineBuffer, InflatedLineBuffer<Cell>>;

/**
 * Line<Cell> API.
 *
 * TODO: Use custom allocator for ensuring cache locality of Cells to sibling lines.
 * TODO: Make the line optimization work.
 */
template <typename Cell>
class Line
{
  public:
    Line() = default;
    Line(Line const&) = default;
    Line(Line&&) noexcept = default;
    Line& operator=(Line const&) = default;
    Line& operator=(Line&&) noexcept = default;

    using TrivialBuffer = TrivialLineBuffer;
    using InflatedBuffer = InflatedLineBuffer<Cell>;
    using Storage = LineStorage<Cell>;
    using value_type = Cell;
    using iterator = typename InflatedBuffer::iterator;
    using reverse_iterator = typename InflatedBuffer::reverse_iterator;
    using const_iterator = typename InflatedBuffer::const_iterator;

    Line(LineFlags _flags, ColumnCount _width, GraphicsAttributes _templateSGR):
        storage_ { TrivialBuffer { _width, _templateSGR } }, flags_ { static_cast<unsigned>(_flags) }
    {
    }

    Line(LineFlags _flags, InflatedBuffer _buffer):
        storage_ { std::move(_buffer) }, flags_ { static_cast<unsigned>(_flags) }
    {
    }

    void reset(LineFlags _flags, GraphicsAttributes _attributes) noexcept
    {
        flags_ = static_cast<unsigned>(_flags);
        if (isTrivialBuffer())
            trivialBuffer().reset(_attributes);
        else
            setBuffer(TrivialBuffer { size(), _attributes });
    }

    void fill(LineFlags _flags,
              GraphicsAttributes const& _attributes,
              char32_t _codepoint,
              uint8_t _width) noexcept
    {
        if (_codepoint == 0)
            reset(_flags, _attributes);
        else
        {
            flags_ = static_cast<unsigned>(_flags);
            for (Cell& cell: inflatedBuffer())
            {
                cell.reset();
                cell.write(_attributes, _codepoint, _width);
            }
        }
    }

    /// Tests if all cells are empty.
    [[nodiscard]] bool empty() const noexcept
    {
        if (isTrivialBuffer())
            return trivialBuffer().text.empty();

        for (auto const& cell: inflatedBuffer())
            if (!cell.empty())
                return false;
        return true;
    }

    /**
     * Fills this line with the given content.
     *
     * @p _start offset into this line of the first charater
     * @p _sgr graphics rendition for the line starting at @c _start until the end
     * @p _ascii the US-ASCII characters to fill with
     */
    void fill(ColumnOffset _start, GraphicsAttributes const& _sgr, std::string_view _ascii)
    {
        auto& buffer = inflatedBuffer();

        assert(unbox<size_t>(_start) + _ascii.size() <= buffer.size());

        auto constexpr ASCII_Width = 1;
        auto const* s = _ascii.data();

        Cell* i = &buffer[unbox<size_t>(_start)];
        Cell* e = i + _ascii.size();
        while (i != e)
            (i++)->write(_sgr, static_cast<char32_t>(*s++), ASCII_Width);

        auto const e2 = buffer.data() + buffer.size();
        while (i != e2)
            (i++)->reset();
    }

    [[nodiscard]] ColumnCount size() const noexcept
    {
        if (isTrivialBuffer())
            return trivialBuffer().displayWidth;
        else
            return ColumnCount::cast_from(inflatedBuffer().size());
    }

    void resize(ColumnCount _count);

    gsl::span<Cell const> trim_blank_right() const noexcept;

    gsl::span<Cell const> cells() const noexcept { return inflatedBuffer(); }

    gsl::span<Cell> useRange(ColumnOffset _start, ColumnCount _count) noexcept
    {
#if defined(__clang__) && __clang_major__ <= 11
        auto const bufferSpan = gsl::span(inflatedBuffer());
        return bufferSpan.subspan(unbox<size_t>(_start), unbox<size_t>(_count));
#else
        // Clang <= 11 cannot deal with this (e.g. FreeBSD 13 defaults to Clang 11).
        return gsl::span(inflatedBuffer()).subspan(unbox<size_t>(_start), unbox<size_t>(_count));
#endif
    }

    Cell& useCellAt(ColumnOffset _column) noexcept
    {
        Require(ColumnOffset(0) <= _column);
        Require(_column <= ColumnOffset::cast_from(size())); // Allow off-by-one for sentinel.
        return inflatedBuffer()[unbox<size_t>(_column)];
    }

    [[nodiscard]] uint8_t cellEmptyAt(ColumnOffset column) const noexcept
    {
        if (isTrivialBuffer())
        {
            Require(ColumnOffset(0) <= column);
            Require(column < ColumnOffset::cast_from(size()));
            return unbox<size_t>(column) >= trivialBuffer().text.size()
                   || trivialBuffer().text[column.as<size_t>()] == 0x20;
        }
        return inflatedBuffer().at(unbox<size_t>(column)).empty();
    }

    [[nodiscard]] uint8_t cellWithAt(ColumnOffset column) const noexcept
    {
        if (isTrivialBuffer())
        {
            Require(ColumnOffset(0) <= column);
            Require(column < ColumnOffset::cast_from(size()));
            return 1; // TODO: When trivial line is to support Unicode, this should be adapted here.
        }
        return inflatedBuffer().at(unbox<size_t>(column)).width();
    }

    [[nodiscard]] LineFlags flags() const noexcept
    {
        return static_cast<LineFlags>(flags_);
    }

    [[nodiscard]] bool marked() const noexcept
    {
        return isFlagEnabled(LineFlags::Marked);
    }
    void setMarked(bool _enable)
    {
        setFlag(LineFlags::Marked, _enable);
    }

    [[nodiscard]] bool wrapped() const noexcept
    {
        return isFlagEnabled(LineFlags::Wrapped);
    }
    void setWrapped(bool _enable)
    {
        setFlag(LineFlags::Wrapped, _enable);
    }

    [[nodiscard]] bool wrappable() const noexcept
    {
        return isFlagEnabled(LineFlags::Wrappable);
    }
    void setWrappable(bool _enable)
    {
        setFlag(LineFlags::Wrappable, _enable);
    }

    [[nodiscard]] LineFlags wrappableFlag() const noexcept
    {
        return wrappable() ? LineFlags::Wrappable : LineFlags::None;
    }
    [[nodiscard]] LineFlags wrappedFlag() const noexcept
    {
        return marked() ? LineFlags::Wrapped : LineFlags::None;
    }
    [[nodiscard]] LineFlags markedFlag() const noexcept
    {
        return marked() ? LineFlags::Marked : LineFlags::None;
    }

    [[nodiscard]] LineFlags inheritableFlags() const noexcept
    {
        auto constexpr Inheritables = unsigned(LineFlags::Wrappable) | unsigned(LineFlags::Marked);
        return static_cast<LineFlags>(flags_ & Inheritables);
    }

    void setFlag(LineFlags _flag, bool _enable) noexcept
    {
        if (_enable)
            flags_ |= static_cast<unsigned>(_flag);
        else
            flags_ &= ~static_cast<unsigned>(_flag);
    }

    [[nodiscard]] bool isFlagEnabled(LineFlags _flag) const noexcept
    {
        return (flags_ & static_cast<unsigned>(_flag)) != 0;
    }

    [[nodiscard]] InflatedBuffer reflow(ColumnCount _newColumnCount);
    [[nodiscard]] std::string toUtf8() const;
    [[nodiscard]] std::string toUtf8Trimmed() const;

    // Returns a reference to this mutable grid-line buffer.
    //
    // If this line has been stored in an optimized state, then
    // the line will be first unpacked into a vector of grid cells.
    InflatedBuffer& inflatedBuffer();
    InflatedBuffer const& inflatedBuffer() const;

    [[nodiscard]] TrivialBuffer& trivialBuffer() noexcept
    {
        return std::get<TrivialBuffer>(storage_);
    }
    [[nodiscard]] TrivialBuffer const& trivialBuffer() const noexcept
    {
        return std::get<TrivialBuffer>(storage_);
    }

    [[nodiscard]] bool isTrivialBuffer() const noexcept
    {
        return std::holds_alternative<TrivialBuffer>(storage_);
    }
    [[nodiscard]] bool isInflatedBuffer() const noexcept
    {
        return !std::holds_alternative<TrivialBuffer>(storage_);
    }

    void setBuffer(TrivialBuffer const& buffer) noexcept
    {
        storage_ = buffer;
    }
    void setBuffer(InflatedBuffer buffer)
    {
        storage_ = std::move(buffer);
    }

    void reset(GraphicsAttributes textAttributes,
               GraphicsAttributes fillAttributes,
               HyperlinkId hyperlink,
               crispy::BufferFragment text,
               ColumnCount columnsUsed)
    {
        storage_ =
            TrivialBuffer { size(), textAttributes, fillAttributes, hyperlink, columnsUsed, std::move(text) };
    }

  private:
    Storage storage_;
    unsigned flags_ = 0;
};

constexpr LineFlags operator|(LineFlags a, LineFlags b) noexcept
{
    return LineFlags(unsigned(a) | unsigned(b));
}

constexpr LineFlags operator~(LineFlags a) noexcept
{
    return LineFlags(~unsigned(a));
}

constexpr LineFlags operator&(LineFlags a, LineFlags b) noexcept
{
    return LineFlags(unsigned(a) & unsigned(b));
}

template <typename Cell>
inline typename Line<Cell>::InflatedBuffer& Line<Cell>::inflatedBuffer()
{
    if (std::holds_alternative<TrivialBuffer>(storage_))
        storage_ = inflate<Cell>(std::get<TrivialBuffer>(storage_));
    return std::get<InflatedBuffer>(storage_);
}

template <typename Cell>
inline typename Line<Cell>::InflatedBuffer const& Line<Cell>::inflatedBuffer() const
{
    return const_cast<Line<Cell>*>(this)->inflatedBuffer();
}

} // namespace terminal

namespace fmt // {{{
{
template <>
struct formatter<terminal::LineFlags>
{
    template <typename ParseContext>
    auto parse(ParseContext& ctx)
    {
        return ctx.begin();
    }
    template <typename FormatContext>
    auto format(const terminal::LineFlags _flags, FormatContext& ctx) const
    {
        static const std::array<std::pair<terminal::LineFlags, std::string_view>, 3> nameMap = {
            std::pair { terminal::LineFlags::Wrappable, std::string_view("Wrappable") },
            std::pair { terminal::LineFlags::Wrapped, std::string_view("Wrapped") },
            std::pair { terminal::LineFlags::Marked, std::string_view("Marked") },
        };
        std::string s;
        for (auto const& mapping: nameMap)
        {
            if ((mapping.first & _flags) != terminal::LineFlags::None)
            {
                if (!s.empty())
                    s += ",";
                s += mapping.second;
            }
        }
        return fmt::format_to(ctx.out(), "{}", s);
    }
};
} // namespace fmt

max_stars_repo_name: christianparpart/libterminal

max_issues_repo_name: christianparpart/libterminal

max_forks_repo_name: christianparpart/libterminal


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

#include <omp.h>

#include "util.h"

int main(int argc, char *argv[]) {
    char* input_fileName1 = argv[1];
    char* input_fileName2 = argv[2];
    int N_doc_bg = atoi(argv[3]);
    int N_kw = atoi(argv[4]);
    int N_obs = atoi(argv[5]);
    int N_iter = atoi(argv[6]);
    double lambda = atof(argv[7]);
    char* output_fileName = argv[8];
    int N_doc = 480000/2;

    int* matrix = (int*) malloc(sizeof(int) * N_obs * N_obs);
    int* matrix_bgi = (int*) malloc(sizeof(int) * N_kw * N_kw);
    int* matrix_padded = (int*) malloc(sizeof(int) * N_obs * N_obs);
    int* true_index = (int*) malloc(sizeof(int) * N_kw);
    int* permutation = (int*) malloc(sizeof(int) * N_obs);
    gsl_matrix* matrix_obs;

    for (int round = 0; round < 10; round++)
    {
        char input_fileName1_extend[40];
        char input_fileName2_extend[40];
        sprintf(input_fileName1_extend, "%s%d", input_fileName1, round);
        sprintf(input_fileName2_extend, "%s%d", input_fileName2, round);

        // Setup
        struct timeval tv1,tv2;
        gettimeofday(&tv1, NULL);
        read_matrix(&true_index, &matrix_bgi, 1.0*N_doc/N_doc_bg, N_kw, input_fileName2_extend);
        read_matrix(&true_index, &matrix, 1.0, N_obs, input_fileName1_extend);
        gettimeofday(&tv2, NULL);
        printf("Reading done: %f.\n", (double) (tv2.tv_usec - tv1.tv_usec) / 1000000 + (double) (tv2.tv_sec - tv1.tv_sec));
        fflush(stdout);

        for (int iter = 0; iter < 10; iter++)
        {
            printf("Run %d\n", iter);
            matrix_obs = gsl_matrix_alloc(N_obs, N_obs);

            gettimeofday(&tv1, NULL);
            pad_matrix(&matrix_padded, &matrix, lambda, N_obs, N_doc);
            gettimeofday(&tv2, NULL);
            printf("Padding done: %f.\n", (double) (tv2.tv_usec - tv1.tv_usec) / 1000000 + (double) (tv2.tv_sec - tv1.tv_sec));
            fflush(stdout);

            gettimeofday(&tv1, NULL);
            observe_matrix(matrix_obs, &matrix_padded, N_obs);
            gettimeofday(&tv2, NULL);
            printf("Observed matrix generated: %f.\n", (double) (tv2.tv_usec - tv1.tv_usec) / 1000000 + (double) (tv2.tv_sec - tv1.tv_sec));
            fflush(stdout);

            gettimeofday(&tv1, NULL);
            attack(matrix_obs, &matrix_bgi, &permutation, N_kw, N_obs, N_doc, N_iter);
            gettimeofday(&tv2, NULL);
            printf("Main attack done: %d.\n",  (tv2.tv_sec - tv1.tv_sec));
            fflush(stdout);

            char output_fileName_full[40];
            sprintf(output_fileName_full, "%s%d-%d", output_fileName, round, iter);
            print_result(output_fileName_full, &permutation, &true_index, N_obs);
            //sprintf(output_fileName_full, "%s%d-%d-full", output_fileName, round, iter);
            //print_full_result(output_fileName_full, &permutation, &true_index, N_obs);
        }
    }

    free(matrix);
    free(matrix_padded);
    gsl_matrix_free(matrix_obs);
    return(0);
}

double log_score(int idx1, int idx2, gsl_matrix* matrix_obs, int** matrix, int** permutation, int N_kw, int N_doc)
{
    if (idx1 == idx2)
        return(0.0);

    int idx1_m = (*permutation)[idx1];
    int idx2_m = (*permutation)[idx2];

    double mean1, mean2, mean3;
    mean3 = 0;
    double var = 0;

    if (gsl_matrix_get(matrix_obs, idx1, idx1) > (*matrix)[idx1_m*N_kw + idx1_m])
    {
        mean1 = (*matrix)[idx1_m*N_kw + idx2_m];
        double x1, x2;
        x1 = gsl_matrix_get(matrix_obs, idx1, idx1) - (*matrix)[idx1_m*N_kw + idx1_m];
        x2 = gsl_matrix_get(matrix_obs, idx2, idx2) - (*matrix)[idx1_m*N_kw + idx2_m];
        mean3 = x1 / N_doc * x2;
        var += x1 / N_doc * x2 * (N_doc - x2) / N_doc;
    }
    else
    {
        mean1 = (*matrix)[idx1_m*N_kw + idx2_m] * gsl_matrix_get(matrix_obs, idx1, idx1) / (*matrix)[idx1_m*N_kw + idx1_m];
        var   += mean1 / (*matrix)[idx1_m*N_kw + idx1_m] * ((*matrix)[idx1_m*N_kw + idx1_m] - gsl_matrix_get(matrix_obs, idx1, idx1));
    }

    if (gsl_matrix_get(matrix_obs, idx2, idx2) > (*matrix)[idx2_m*N_kw + idx2_m])
    {
        mean2 = mean1;
        double x1, x2;
        x1 = gsl_matrix_get(matrix_obs, idx2, idx2) - (*matrix)[idx2_m*N_kw + idx2_m];
        x2 = gsl_matrix_get(matrix_obs, idx1, idx1) - (*matrix)[idx1_m*N_kw + idx2_m];
        mean3 += x1 / N_doc * x2;
        var += x1 / N_doc * x2 * (N_doc - x2) / N_doc;
    }
    else
    {
        mean2 = mean1 * gsl_matrix_get(matrix_obs, idx2, idx2) / (*matrix)[idx2_m*N_kw + idx2_m];
        var  += mean2 / (*matrix)[idx2_m*N_kw + idx2_m] * ((*matrix)[idx2_m*N_kw + idx2_m] - gsl_matrix_get(matrix_obs, idx2, idx2));
    }

    var += 1.0 * (*matrix)[idx1_m*N_kw + idx2_m] / N_doc * (N_doc - (*matrix)[idx2_m*N_kw + idx2_m]);

    double score = gsl_ran_gaussian_pdf(mean2 + mean3 - gsl_matrix_get(matrix_obs, idx1, idx2), sqrt(var));

    if (score == 0)
        return(-500.0);
    return(log(score));
}

void attack(gsl_matrix* matrix_obs, int** matrix, int** permutation, int N_kw, int N_obs, int N_doc, int N_iter)
{
    // Initialise data structures
    double* score_matrix = (double*) malloc(sizeof(double) * N_obs * N_obs);
    double* score_row1 = (double*) malloc(sizeof(double) * N_obs);
    double* score_row2 = (double*) malloc(sizeof(double) * N_obs);

    int* permutation_tmp = (int*) malloc(sizeof(int) * N_obs);
    int* permutation_inv = (int*) malloc(sizeof(int) * N_kw);

    // Initialise permutations
    for (int ii = 0; ii < N_obs; ii++)
        (*permutation)[ii] = ii;
    for (int ii = 0; ii < N_obs; ii++)
        permutation_tmp[ii] = ii;
    for (int ii = 0; ii < N_obs; ii++)
        permutation_inv[ii] = ii;
    for (int ii = N_obs; ii < N_kw; ii++)
        permutation_inv[ii] = -1;

    // Initialising RNG
    const gsl_rng_type * T;
    gsl_rng * r;
    gsl_rng_env_setup();
    T = gsl_rng_default;
    r = gsl_rng_alloc (T);

    struct timeval tv1,tv2;
    gettimeofday(&tv1, NULL);

    // Compute initial score
    #pragma omp parallel for shared(score_matrix, matrix_obs, matrix)
    for (int ii = 0; ii < N_obs * N_obs; ii++)
        score_matrix[ii] = log_score((int) (ii / N_obs), ii % N_obs, matrix_obs, matrix, permutation, N_kw, N_doc);
    gettimeofday(&tv2, NULL);
    printf("Initial score computed: %f.\n", (double) (tv2.tv_usec - tv1.tv_usec) / 1000000 + (double) (tv2.tv_sec - tv1.tv_sec));

    // Iterations of simulated annealing
    double temp = (double) N_kw;
    int N_stuck = 0;
    for (int iter = 0; iter < N_iter; iter++)
    {
	    /* Status code */
	    if (iter % (N_iter / 10) == 0)
        {
            gettimeofday(&tv1, NULL);
            printf("Iteration: %d, %d, %d.\n", iter, N_stuck, (int) (tv1.tv_sec - tv2.tv_sec));
            fflush(stdout);
            gettimeofday(&tv2, NULL);
        }

        if (N_stuck >= 20000)
            iter = N_iter;

	    /* Main code */
        int idx1, idx2;
        permutation_generation(&idx1, &idx2, &permutation_tmp, permutation, &permutation_inv, N_kw, N_obs);

        int ii = 0;
        #pragma omp parallel for shared(score_row1)
        for (ii = 0; ii < N_obs; ii++)
            score_row1[ii] = log_score(idx1, ii, matrix_obs, matrix, &permutation_tmp, N_kw, N_doc);

        if (idx2 >= 0)
            #pragma omp parallel for shared(score_row2)
            for (ii = 0; ii < N_obs; ii++)
                score_row2[ii] = log_score(idx2, ii, matrix_obs, matrix, &permutation_tmp, N_kw, N_doc);

        double score_diff = 0;
        for (int ii = 0; ii < N_obs; ii++)
            score_diff += score_row1[ii];
        for (int ii = 0; ii < N_obs; ii++)
            score_diff -= score_matrix[idx1*N_obs + ii];
        if (idx2 >= 0)
        {
            for (int ii = 0; ii < N_obs; ii++)
                score_diff += score_row2[ii];
            for (int ii = 0; ii < N_obs; ii++)
                score_diff -= score_matrix[idx2*N_obs + ii];
        }

        // compute difference in score, with exponentiation
        score_diff = score_diff / temp;

        if (score_diff < -40)
            score_diff = 0;
        else if (score_diff > 0)
            score_diff = 1.01;
        else
            score_diff = exp(score_diff);

        if (gsl_ran_flat(r, 0, 1) < score_diff)
        {
            // Update the scores
            for (int ii = 0; ii < N_obs; ii++)
                score_matrix[idx1*N_obs + ii] = score_row1[ii];
            for (int ii = 0; ii < N_obs; ii++)
                score_matrix[ii*N_obs + idx1] = score_row1[ii];
            if (idx2 >= 0)
            {
                for (int ii = 0; ii < N_obs; ii++)
                    score_matrix[idx2*N_obs + ii] = score_row2[ii];
                for (int ii = 0; ii < N_obs; ii++)
                    score_matrix[ii*N_obs + idx2] = score_row2[ii];
            }

            // Update the permutation
            permutation_inv[(*permutation)[idx1]] = -1;
            (*permutation)[idx1] = permutation_tmp[idx1];
            permutation_inv[permutation_tmp[idx1]] = idx1;
            if (idx2 >= 0)
            {
                (*permutation)[idx2] = permutation_tmp[idx2];
                permutation_inv[permutation_tmp[idx2]] = idx2;
            }
            N_stuck = 0;
        }
        else
        {
            // Update the permutation
            permutation_tmp[idx1] = (*permutation)[idx1];
            if (idx2 >= 0)
                permutation_tmp[idx2] = (*permutation)[idx2];
            N_stuck += 1;
        }

        temp *= 0.995;
    }

    free(score_matrix);
    free(score_row1);
    free(score_row2);
    gsl_rng_free(r);
}

void print_result(char* output_fileName, int** permutation, int** true_index, int N_obs)
{
    FILE* fp = fopen(output_fileName, "w");
    int count = 0;
    for (int ii = 0; ii < N_obs; ii++)
        if ((*permutation)[ii] == (*true_index)[ii])
            count++;

    fprintf(fp, "%d\n", count);
    fclose(fp);
    printf("Success: %d/%d.\n", count, N_obs);
}

void print_full_result(char* output_fileName, int** permutation, int** true_index, int N_obs)
{
    FILE* fp = fopen(output_fileName, "w");
    int count = 0;
    for (int ii = 0; ii < N_obs; ii++)
        fprintf(fp, "%d, %d\n", (*permutation)[ii], (*true_index)[ii]);
    fclose(fp);
}

ext: c

max_stars_repo_name: RethinkingSSE/Attacks-on-SSE

max_issues_repo_name: RethinkingSSE/Attacks-on-SSE

max_forks_repo_name: RethinkingSSE/Attacks-on-SSE


text:
/* multifit_nlinear/subspace2D.c
 *
 * Copyright (C) 2016 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 <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_multifit_nlinear.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_poly.h>

/*
 * This module implements a 2D subspace trust region subproblem method,
 * as outlined in
 *
 * [1] G. A. Shultz, R. B. Schnabel, and R. H. Byrd
 *     A Family of Trust-Region-Based Algorithms for Unconstrained
 *     Minimization with Strong Global Convergence Properties,
 *     SIAM Journal on Numerical Analysis 1985 22:1, 47-67
 *
 * [2] R. H. Byrd, R. B. Schnabel, G. A. Shultz,
 *     Approximate solution of the trust region problem by
 *     minimization over two-dimensional subspaces,
 *     Mathematical Programming, January 1988, Volume 40,
 *     Issue 1, pp 247-263
 *
 * The idea is to solve:
 *
 * min_{dx} g^T dx + 1/2 dx^T B dx
 *
 * with constraints:
 *
 * ||D dx|| <= delta
 * dx \in span{dx_sd, dx_gn}
 *
 * where B is the Hessian matrix, B = J^T J
 *
 * The steps are as follows:
 *
 * 1. preloop:
 *    a. Compute Gauss-Newton and steepest descent vectors,
 *       dx_gn, dx_sd
 *    b. Compute an orthonormal basis for span(D dx_sd, D dx_gn) by
 *       constructing W = [ D dx_sd, D dx_gn ] and performing a QR
 *       decomposition of W. The 2 columns of the Q matrix
 *       will then span the column space of W. W should have rank 2
 *       unless D*dx_sd and D*dx_gn are parallel, in which case it will
 *       have rank 1.
 *    c. Precompute various quantities needed for the step calculation
 *
 * 2. step:
 *    a. If the Gauss-Newton step is inside the trust region, use it
 *    b. if W has rank 1, we cannot form a 2D subspace, so in this case
 *       follow the steepest descent direction to the trust region boundary
 *       and use that as the step.
 *    c. In the full rank 2 case, if the GN point is outside the trust region,
 *       then the minimizer of the objective function lies on the trust
 *       region boundary. Therefore the minimization problem becomes:
 *
 *       min_{dx} g^T dx + 1/2 dx^T B dx, with ||dx|| = delta, dx = Q * x
 *
 *       where x is a 2-vector to be determined and the columns of Q are
 *       the orthonormal basis vectors of the subspace. Note the equality
 *       constraint now instead of <=. In terms of the new variable x,
 *       the minimization problem becomes:
 *
 *       min_x subg^T x + 1/2 x^T subB x, with ||Q*x|| = ||x|| = delta
 *
 *       where:
 *         subg = Q^T g   (2-by-1)
 *         subB = Q^T B Q (2-by-2)
 *
 *       This equality constrained 2D minimization problem can be solved
 *       with a Lagrangian multiplier, which results in a 4th degree polynomial
 *       equation to be solved. The equation is:
 *
 *         lambda^4  1
 *       + lambda^3  2 tr(B)
 *       + lambda^2  (tr(B)^2 + 2 det(B) - g^T g / delta^2)
 *       + lambda^1  (2 det(B) tr(B) - 2 g^T adj(B)^T g / delta^2)
 *       + lambda^0  (det(B)^2 - g^T adj(B)^T adj(B) g / delta^2)
 *
 *       where adj(B) is the adjugate matrix of B.
 *
 *       We then check each of the 4 solutions for lambda to determine which
 *       lambda results in the smallest objective function value. This x
 *       is then used to construct the final step: dx = Q*x
 */

typedef struct
{
  size_t n;                  /* number of observations */
  size_t p;                  /* number of parameters */
  gsl_vector *dx_gn;         /* Gauss-Newton step, size p */
  gsl_vector *dx_sd;         /* steepest descent step, size p */
  double norm_Dgn;           /* || D dx_gn || */
  double norm_Dsd;           /* || D dx_sd || */
  gsl_vector *workp;         /* workspace, length p */
  gsl_vector *workn;         /* workspace, length n */
  gsl_matrix *W;             /* orthonormal basis for 2D subspace, p-by-2 */
  gsl_matrix *JQ;            /* J * Q, n-by-p */
  gsl_vector *tau;           /* Householder scalars */
  gsl_vector *subg;          /* subspace gradient = W^T g, 2-by-1 */
  gsl_matrix *subB;          /* subspace Hessian = W^T B W, 2-by-2 */
  gsl_permutation *perm;     /* permutation matrix */

  double trB;                /* Tr(subB) */
  double detB;               /* det(subB) */
  double normg;              /* || subg || */
  double term0;              /* g^T adj(B)^T adj(B) g */
  double term1;              /* g^T adj(B)^T g */

  size_t rank;               /* rank of [ dx_sd, dx_gn ] matrix */

  gsl_poly_complex_workspace *poly_p;

  /* tunable parameters */
  gsl_multifit_nlinear_parameters params;
} subspace2D_state_t;

#include "common.c"

static void * subspace2D_alloc (const void * params, const size_t n, const size_t p);
static void subspace2D_free(void *vstate);
static int subspace2D_init(const void *vtrust_state, void *vstate);
static int subspace2D_preloop(const void * vtrust_state, void * vstate);
static int subspace2D_step(const void * vtrust_state, const double delta,
                           gsl_vector * dx, void * vstate);
static int subspace2D_preduction(const void * vtrust_state, const gsl_vector * dx,
                                 double * pred, void * vstate);
static int subspace2D_solution(const double lambda, gsl_vector * x,
                               subspace2D_state_t * state);
static double subspace2D_objective(const gsl_vector * x, subspace2D_state_t * state);
static int subspace2D_calc_gn(const gsl_multifit_nlinear_trust_state * trust_state, gsl_vector * dx);
static int subspace2D_calc_sd(const gsl_multifit_nlinear_trust_state * trust_state, gsl_vector * dx,
                              subspace2D_state_t * state);

static void *
subspace2D_alloc (const void * params, const size_t n, const size_t p)
{
  const gsl_multifit_nlinear_parameters *par = (const gsl_multifit_nlinear_parameters *) params;
  subspace2D_state_t *state;

  state = calloc(1, sizeof(subspace2D_state_t));
  if (state == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate subspace2D state", GSL_ENOMEM);
    }

  state->dx_gn = gsl_vector_alloc(p);
  if (state->dx_gn == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for dx_gn", GSL_ENOMEM);
    }

  state->dx_sd = gsl_vector_alloc(p);
  if (state->dx_sd == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for dx_sd", GSL_ENOMEM);
    }

  state->workp = gsl_vector_alloc(p);
  if (state->workp == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for workp", GSL_ENOMEM);
    }

  state->workn = gsl_vector_alloc(n);
  if (state->workn == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for workn", GSL_ENOMEM);
    }

  state->W = gsl_matrix_alloc(p, 2);
  if (state->W == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for W", GSL_ENOMEM);
    }

  state->JQ = gsl_matrix_alloc(n, p);
  if (state->JQ == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for JQ", GSL_ENOMEM);
    }

  state->tau = gsl_vector_alloc(2);
  if (state->tau == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for tau", GSL_ENOMEM);
    }

  state->subg = gsl_vector_alloc(2);
  if (state->subg == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for subg", GSL_ENOMEM);
    }

  state->subB = gsl_matrix_alloc(2, 2);
  if (state->subB == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for subB", GSL_ENOMEM);
    }

  state->perm = gsl_permutation_alloc(2);
  if (state->perm == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for perm", GSL_ENOMEM);
    }

  state->poly_p = gsl_poly_complex_workspace_alloc(5);
  if (state->poly_p == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for poly workspace", GSL_ENOMEM);
    }

  state->n = n;
  state->p = p;
  state->rank = 0;
  state->params = *par;

  return state;
}

static void
subspace2D_free(void *vstate)
{
  subspace2D_state_t *state = (subspace2D_state_t *) vstate;

  if (state->dx_gn)
    gsl_vector_free(state->dx_gn);

  if (state->dx_sd)
    gsl_vector_free(state->dx_sd);

  if (state->workp)
    gsl_vector_free(state->workp);

  if (state->workn)
    gsl_vector_free(state->workn);

  if (state->W)
    gsl_matrix_free(state->W);

  if (state->JQ)
    gsl_matrix_free(state->JQ);

  if (state->tau)
    gsl_vector_free(state->tau);

  if (state->subg)
    gsl_vector_free(state->subg);

  if (state->subB)
    gsl_matrix_free(state->subB);

  if (state->perm)
    gsl_permutation_free(state->perm);

  if (state->poly_p)
    gsl_poly_complex_workspace_free(state->poly_p);

  free(state);
}

/*
subspace2D_init()
  Initialize subspace2D solver

Inputs: vtrust_state - trust state
        vstate       - workspace

Return: success/error
*/

static int
subspace2D_init(const void *vtrust_state, void *vstate)
{
  (void)vtrust_state;
  (void)vstate;

  return GSL_SUCCESS;
}

/*
subspace2D_preloop()
  Initialize subspace2D method prior to iteration loop.
This involves computing the Gauss-Newton step and
steepest descent step

Notes: on output,
1) state->dx_gn contains Gauss-Newton step
2) state->dx_sd contains steepest descent step
3) state->rank contains the rank([dx_sd, dx_gn])
4) if full rank subspace (rank = 2), then:
   state->trB = Tr(subB)
   state->detB = det(subB)
   state->normg = || subg ||
*/

static int
subspace2D_preloop(const void * vtrust_state, void * vstate)
{
  int status;
  const gsl_multifit_nlinear_trust_state *trust_state =
    (const gsl_multifit_nlinear_trust_state *) vtrust_state;
  subspace2D_state_t *state = (subspace2D_state_t *) vstate;
  gsl_vector_view v;
  double work_data[2];
  gsl_vector_view work = gsl_vector_view_array(work_data, 2);
  int signum;

  /* calculate Gauss-Newton step */
  status = subspace2D_calc_gn(trust_state, state->dx_gn);
  if (status)
    return status;

  /* now calculate the steepest descent step */
  status = subspace2D_calc_sd(trust_state, state->dx_sd, state);
  if (status)
    return status;

  /* store norms */
  state->norm_Dgn = scaled_enorm(trust_state->diag, state->dx_gn);
  state->norm_Dsd = scaled_enorm(trust_state->diag, state->dx_sd);

  /*
   * now compute orthonormal basis for span(D dx_sd, D dx_gn) using
   * QR decomposition; set W = [ D dx_sd, D dx_gn ] and normalize each
   * column to unit magnitude. Then the Q matrix will form a basis for Col(W)
   */

  v = gsl_matrix_column(state->W, 0);
  gsl_vector_memcpy(&v.vector, state->dx_sd);
  gsl_vector_mul(&v.vector, trust_state->diag);
  if (state->norm_Dsd != 0)
    gsl_vector_scale(&v.vector, 1.0 / state->norm_Dsd);

  v = gsl_matrix_column(state->W, 1);
  gsl_vector_memcpy(&v.vector, state->dx_gn);
  gsl_vector_mul(&v.vector, trust_state->diag);
  if (state->norm_Dgn != 0)
    gsl_vector_scale(&v.vector, 1.0 / state->norm_Dgn);

  /* use a rank revealing QR decomposition in case dx_sd and dx_gn
   * are parallel */
  gsl_linalg_QRPT_decomp(state->W, state->tau, state->perm, &signum, &work.vector);

  /* check for parallel dx_sd, dx_gn, in which case rank will be 1 */
  state->rank = gsl_linalg_QRPT_rank(state->W, -1.0);

  if (state->rank == 2)
    {
      /*
       * full rank subspace, compute:
       * subg = Q^T D^{-1} g
       * subB = Q^T D^{-1} B D^{-1} Q where B = J^T J
       */
      const size_t p = state->p;
      size_t i;
      gsl_matrix_view JQ = gsl_matrix_submatrix(state->JQ, 0, 0, state->n, GSL_MIN(2, p));
      double B00, B10, B11, g0, g1;

      /* compute subg */
      gsl_vector_memcpy(state->workp, trust_state->g);
      gsl_vector_div(state->workp, trust_state->diag);
      gsl_linalg_QR_QTvec(state->W, state->tau, state->workp);

      g0 = gsl_vector_get(state->workp, 0);
      g1 = gsl_vector_get(state->workp, 1);

      gsl_vector_set(state->subg, 0, g0);
      gsl_vector_set(state->subg, 1, g1);

      /* compute subB */

      /* compute J D^{-1} */
      gsl_matrix_memcpy(state->JQ, trust_state->J);

      for (i = 0; i < p; ++i)
        {
          gsl_vector_view c = gsl_matrix_column(state->JQ, i);
          double di = gsl_vector_get(trust_state->diag, i);
          gsl_vector_scale(&c.vector, 1.0 / di);
        }

      /* compute J D^{-1} Q */
      gsl_linalg_QR_matQ(state->W, state->tau, state->JQ);

      /* compute subB = Q^T D^{-1} J^T J D^{-1} Q */
      gsl_blas_dsyrk(CblasLower, CblasTrans, 1.0, &JQ.matrix, 0.0, state->subB);

      B00 = gsl_matrix_get(state->subB, 0, 0);
      B10 = gsl_matrix_get(state->subB, 1, 0);
      B11 = gsl_matrix_get(state->subB, 1, 1);

      state->trB = B00 + B11;
      state->detB = B00*B11 - B10*B10;
      state->normg = gsl_blas_dnrm2(state->subg);

      /* g^T adj(B)^T adj(B) g */
      state->term0 = (B10*B10 + B11*B11)*g0*g0 -
                     2*B10*(B00 + B11)*g0*g1 +
                     (B00*B00 + B10*B10)*g1*g1;

      /* g^T adj(B)^T g */
      state->term1 = B11 * g0 * g0 + g1 * (B00*g1 - 2*B10*g0);

    }

  return GSL_SUCCESS;
}

/*
subspace2D_step()
  Calculate a new step with 2D subspace method. Based on [1]. We
seek a vector dx in span{dx_gn, dx_sd} which minimizes the model
function subject to ||dx|| <= delta
*/

static int
subspace2D_step(const void * vtrust_state, const double delta,
                gsl_vector * dx, void * vstate)
{
  const gsl_multifit_nlinear_trust_state *trust_state =
    (const gsl_multifit_nlinear_trust_state *) vtrust_state;
  subspace2D_state_t *state = (subspace2D_state_t *) vstate;

  if (state->norm_Dgn <= delta)
    {
      /* Gauss-Newton step is inside trust region, use it as final step
       * since it is the global minimizer of the quadratic model function */
      gsl_vector_memcpy(dx, state->dx_gn);
    }
  else if (state->rank < 2)
    {
      /* rank of [dx_sd, dx_gn] is 1, meaning dx_sd and dx_gn
       * are parallel so we can't form a 2D subspace. Follow the steepest
       * descent direction to the trust region boundary as our step */
      gsl_vector_memcpy(dx, state->dx_sd);
      gsl_vector_scale(dx, delta / state->norm_Dsd);
    }
  else
    {
      int status;
      const double delta_sq = delta * delta;
      double u = state->normg / delta;
      double a[5];
      double z[8];

#if 1
      a[0] = state->detB * state->detB - state->term0 / delta_sq;
      a[1] = 2 * state->detB * state->trB - 2 * state->term1 / delta_sq;
      a[2] = state->trB * state->trB + 2 * state->detB - u * u;
      a[3] = 2 * state->trB;
      a[4] = 1.0;
#else
      double TrB_D = state->trB * delta;
      double detB_D = state->detB * delta;
      double normg_sq = state->normg * state->normg;

      a[0] = detB_D * detB_D - state->term0;
      a[1] = 2 * state->detB * state->trB * delta_sq - 2 * state->term1;
      a[2] = TrB_D * TrB_D + 2 * state->detB * delta_sq - normg_sq;
      a[3] = 2 * state->trB * delta_sq;
      a[4] = delta_sq;
#endif

      status = gsl_poly_complex_solve(a, 5, state->poly_p, z);
      if (status == GSL_SUCCESS)
        {
          size_t i;
          double min = 0.0;
          int mini = -1;
          double x_data[2];
          gsl_vector_view x = gsl_vector_view_array(x_data, 2);

          /*
           * loop through all four values of the Lagrange multiplier
           * lambda. For each lambda, evaluate the objective function
           * with Re(lambda) to determine which lambda minimizes the
           * function
           */
          for (i = 0; i < 4; ++i)
            {
              double cost, normx;

              /*fprintf(stderr, "root: %.12e + %.12e i\n",
                      z[2*i], z[2*i+1]);*/

              status = subspace2D_solution(z[2*i], &x.vector, state);
              if (status != GSL_SUCCESS)
                continue; /* singular matrix system */

              /* ensure ||x|| = delta */

              normx = gsl_blas_dnrm2(&x.vector);
              if (normx == 0.0)
                continue;

              gsl_vector_scale(&x.vector, delta / normx);

              /* evaluate objective function to determine minimizer */
              cost = subspace2D_objective(&x.vector, state);
              if (mini < 0 || cost < min)
                {
                  mini = (int) i;
                  min = cost;
                }
            }

          if (mini < 0)
            {
              /* did not find minimizer - should not get here */
              return GSL_FAILURE;
            }
          else
            {
              /* compute x which minimizes objective function */
              subspace2D_solution(z[2*mini], &x.vector, state);

              /* dx = Q * x */
              gsl_vector_set_zero(dx);
              gsl_vector_set(dx, 0, gsl_vector_get(&x.vector, 0));
              gsl_vector_set(dx, 1, gsl_vector_get(&x.vector, 1));
              gsl_linalg_QR_Qvec(state->W, state->tau, dx);

              /* compute final dx by multiplying by D^{-1} */
              gsl_vector_div(dx, trust_state->diag);
            }
        }
      else
        {
          GSL_ERROR ("gsl_poly_complex_solve failed", status);
        }
    }

  return GSL_SUCCESS;
}

static int
subspace2D_preduction(const void * vtrust_state, const gsl_vector * dx,
                  double * pred, void * vstate)
{
  const gsl_multifit_nlinear_trust_state *trust_state =
    (const gsl_multifit_nlinear_trust_state *) vtrust_state;
  subspace2D_state_t *state = (subspace2D_state_t *) vstate;

  *pred = quadratic_preduction(trust_state->f, trust_state->J, dx, state->workn);

  return GSL_SUCCESS;
}

/* solve 2D subspace problem: (B + lambda*I) x = -g */
static int
subspace2D_solution(const double lambda, gsl_vector * x,
                    subspace2D_state_t * state)
{
  int status = GSL_SUCCESS;
  double C_data[4];
  gsl_matrix_view C = gsl_matrix_view_array(C_data, 2, 2);
  double B00 = gsl_matrix_get(state->subB, 0, 0);
  double B10 = gsl_matrix_get(state->subB, 1, 0);
  double B11 = gsl_matrix_get(state->subB, 1, 1);

  /* construct C = B + lambda*I */
  gsl_matrix_set(&C.matrix, 0, 0, B00 + lambda);
  gsl_matrix_set(&C.matrix, 1, 0, B10);
  gsl_matrix_set(&C.matrix, 0, 1, B10);
  gsl_matrix_set(&C.matrix, 1, 1, B11 + lambda);

  /* use modified Cholesky in case C is not positive definite */
  gsl_linalg_mcholesky_decomp(&C.matrix, state->perm, NULL);
  gsl_linalg_mcholesky_solve(&C.matrix, state->perm, state->subg, x);

  gsl_vector_scale(x, -1.0);

  return status;
}

/* evaluate 2D objective function: f(x) = g^T x + 1/2 x^T B x */
static double
subspace2D_objective(const gsl_vector * x, subspace2D_state_t * state)
{
  double f;
  double y_data[2];
  gsl_vector_view y = gsl_vector_view_array(y_data, 2);

  /* compute: y = g + 1/2 B x */
  gsl_vector_memcpy(&y.vector, state->subg);
  gsl_blas_dsymv(CblasLower, 0.5, state->subB, x, 1.0, &y.vector);

  /* compute: f = x^T ( g + 1/2 B x ) */
  gsl_blas_ddot(x, &y.vector, &f);

  return f;
}

/*
subspace2D_calc_gn()
  Calculate Gauss-Newton step which satisfies:

J dx_gn = -f

Inputs: trust_state - trust state variables
        dx          - (output) Gauss-Newton step

Return: success/error
*/

static int
subspace2D_calc_gn(const gsl_multifit_nlinear_trust_state * trust_state, gsl_vector * dx)
{
  int status;
  const gsl_multifit_nlinear_parameters *params = trust_state->params;

  /* initialize linear least squares solver */
  status = (params->solver->init)(trust_state, trust_state->solver_state);
  if (status)
    return status;

  /* prepare the linear solver to compute Gauss-Newton step */
  status = (params->solver->presolve)(0.0, trust_state, trust_state->solver_state);
  if (status)
    return status;

  /* solve: J dx_gn = -f for Gauss-Newton step */
  status = (params->solver->solve)(trust_state->f,
                                   dx,
                                   trust_state,
                                   trust_state->solver_state);
  if (status)
    return status;

  return GSL_SUCCESS;
}

/*
subspace2D_calc_sd()
  Calculate steepest descent step,

dx_sd = - || D^{-1} g ||^2 / || J D^{-2} g ||^2 D^{-2} g

Inputs: trust_state - trust state variables
        dx          - (output) steepest descent vector
        state       - workspace

Return: success/error
*/

static int
subspace2D_calc_sd(const gsl_multifit_nlinear_trust_state * trust_state, gsl_vector * dx,
                   subspace2D_state_t * state)
{
  double norm_Dinvg;   /* || D^{-1} g || */
  double norm_JDinv2g; /* || J D^{-2} g || */
  double alpha;        /* || D^{-1} g ||^2 / || J D^{-2} g ||^2 */
  double u;

  /* compute workp = D^{-1} g and its norm */
  gsl_vector_memcpy(state->workp, trust_state->g);
  gsl_vector_div(state->workp, trust_state->diag);
  norm_Dinvg = gsl_blas_dnrm2(state->workp);

  /* compute workp = D^{-2} g */
  gsl_vector_div(state->workp, trust_state->diag);

  /* compute: workn = J D^{-2} g */
  gsl_blas_dgemv(CblasNoTrans, 1.0, trust_state->J, state->workp, 0.0, state->workn);
  norm_JDinv2g = gsl_blas_dnrm2(state->workn);

  u = norm_Dinvg / norm_JDinv2g;
  alpha = u * u;

  /* dx_sd = -alpha D^{-2} g */
  gsl_vector_memcpy(dx, state->workp);
  gsl_vector_scale(dx, -alpha);

  return GSL_SUCCESS;
}

static const gsl_multifit_nlinear_trs subspace2D_type =
{
  "2D-subspace",
  subspace2D_alloc,
  subspace2D_init,
  subspace2D_preloop,
  subspace2D_step,
  subspace2D_preduction,
  subspace2D_free
};

const gsl_multifit_nlinear_trs *gsl_multifit_nlinear_trs_subspace2D = &subspace2D_type;

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/* Copyright (c) 2014, Giuseppe Argentieri <giuseppe.argentieri@ts.infn.it>

 * All rights reserved.
 *
 * 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 THE COPYRIGHT HOLDERS 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 COPYRIGHT HOLDER 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.
 *
 *
 *
 *       Filename:  evol.c
 *
 *    Description:  Evolution of the density matrix with a given generator
 *
 *        Version:  1.0
 *        Created:  02/05/2014 14:40:25
 *       Revision:  none
 *        License:  BSD
 *
 *         Author:  Giuseppe Argentieri (ga), giuseppe.argentieri@ts.infn.it
 *   Organization:  Università degli Studi di Trieste
 *
 *
 */

#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>

#include "funcs.h"

/*
 *      FUNCTION
 *         Name:  generator
 *  Description:  Setting the dydt in Bloch form
 *
 */
int generator ( double t, const double y[], double dydt[], void* PARS )
{
	/* Extracting the matrix address from PARS */
	gsl_matrix* M = (gsl_matrix*) PARS ;

	/* Generator */
	unsigned int i , j ;
	for ( i = 0; i < 4 ; i++ )
		dydt[i] = 0 ;

	for ( i = 1 ; i < 4 ; i++ )
		for ( j = 0 ; j < 4 ; j++ )
			dydt[i] = dydt[i] + gsl_matrix_get( M, i, j )*y[j] ;

	return GSL_SUCCESS;
}		/* -----  end of function generator  ----- */

/*
 *      FUNCTION
 *         Name:  jac
 *  Description:  Jacobian matrix of the generator (which is exactly the Bloch matrix)
 *
 */
int jac ( double t, const double y[], double dfdy[] , double dfdt[], void* PARS )
{
	/* Taking the generator Bloch's matrix address from PARS */
	gsl_matrix* bloch = (gsl_matrix*) PARS ;

	/*  Creating a matrix view of the array dfdy */
	gsl_matrix_view m = gsl_matrix_view_array ( dfdy , 4, 4 ) ;

	/* Initializing the jacobian matrix with the Bloch generator */
	unsigned int i, j ;
	for ( i = 0 ; i < 4 ; i++ )
		for ( j = 0 ; j < 4 ; j++ )
			gsl_matrix_set ( &m.matrix, i, j, gsl_matrix_get(bloch,i,j) ) ;

	for ( i = 0 ; i < 4 ; i++ )
		dfdt[i] = 0 ;

	return GSL_SUCCESS;
}		/* -----  end of function jac  ----- */

/*
 *      FUNCTION
 *         Name:  evol
 *  Description:
 *
 */
int evol ( double t, gsl_vector* state, double step,
		gsl_odeiv2_evolve* e, gsl_odeiv2_control* c, gsl_odeiv2_step* s,
		gsl_odeiv2_system* sys )
{
	/* Creating the array 'rho' for the vector 'state' */
	double rho[4] ;
	unsigned int i ;
	for ( i = 0 ; i < 4 ; i++ )
		rho[i] = VECTOR( state, i ) ;

	/*  evolving the array 'rho'  */
	int status = gsl_odeiv2_evolve_apply_fixed_step	( e, c, s, sys, &t, step, rho ) ;
	if ( status != GSL_SUCCESS )
	{
		printf("STATUS: %d\n", status) ;
		exit(1) ;
	}

	/* updating vector 'state' */
	for ( i = 0 ; i < 4 ; i++ )
		gsl_vector_set( state, i, rho[i] ) ;

	return status ;
}		/* -----  end of function evol  ----- */


text:
/* -----------------------------------------------------------------------------
 * Copyright 2021 Jonathan Haigh
 * SPDX-License-Identifier: MIT
 * ---------------------------------------------------------------------------*/

#ifndef SQ_INCLUDE_GUARD_core_Token_h_
#define SQ_INCLUDE_GUARD_core_Token_h_

#include "Token.fwd.h"

#include "core/typeutil.h"

#include <gsl/gsl>
#include <iosfwd>
#include <regex>
#include <string_view>
#include <vector>

namespace sq {

enum class TokenKind : int {
  BoolFalse,
  BoolTrue,
  Colon,
  Comma,
  Dot,
  DQString,
  Eof,
  Equals,
  Float,
  GreaterThan,
  GreaterThanOrEqualTo,
  Identifier,
  Integer,
  LBrace,
  LBracket,
  LessThan,
  LessThanOrEqualTo,
  LParen,
  RBrace,
  RBracket,
  RParen
};

class Token {
public:
  /**
   * Create a Token object.
   *
   * @param query the full query string in which the token was found.
   * @param pos the character position within the query at which the token
   *        was found.
   * @param len the length, in characters, of the token.
   * @param kind the kind of the token.
   */
  Token(std::string_view query, gsl::index pos, gsl::index len,
        TokenKind kind) noexcept;

  /**
   * Get the full query string in which the token was found.
   */
  SQ_ND std::string_view query() const noexcept;

  /**
   * Get the character position within the query at which the token was
   * found.
   */
  SQ_ND gsl::index pos() const noexcept;

  /**
   * Get the length, in characters, of the token.
   */
  SQ_ND gsl::index len() const noexcept;

  /**
   * Get a std::string_view pointing to the characters of the token.
   */
  SQ_ND std::string_view view() const noexcept;

  /**
   * Get the kind of the token.
   */
  SQ_ND TokenKind kind() const noexcept;

private:
  std::string_view query_;
  gsl::index pos_;
  gsl::index len_;
  TokenKind kind_;
};

std::ostream &operator<<(std::ostream &os, TokenKind kind);

/**
 * Print information about a token.
 *
 * The text printed includes information about the position of the token in the
 * input query, not just the characters that make up the token.
 */
std::ostream &operator<<(std::ostream &os, const Token &token);

} // namespace sq

#endif // SQ_INCLUDE_GUARD_core_Token_h_

ext: h

max_stars_repo_name: jonathanhaigh/sq

max_issues_repo_name: jonathanhaigh/sq

max_forks_repo_name: jonathanhaigh/sq


text:
// Copyright 2019 John McFarlane
//
// 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 WSS_BOARD_H
#define WSS_BOARD_H

#include "coord.h"
#include "grid.h"
#include "load_buffer.h"

#include <ssize.h>
#include <wss_assert.h>

#include <fmt/printf.h>
#include <gsl/string_span>

#include <array>
#include <optional>
#include <unordered_map>

template<typename T>
class board {
public:
    board(board&&) = default;

    explicit board(int init_edge)
            :edge{init_edge}, cells{std::make_unique<T[]>(edge*edge)} { }

    int size() const
    {
        return edge;
    }

    T const& cell(coord c) const
    {
        WSS_ASSERT(c[0]>=0);
        WSS_ASSERT(c[0]<edge);
        WSS_ASSERT(c[1]>=0);
        WSS_ASSERT(c[1]<edge);

        return cells.get()[c[0]+c[1]*edge];
    }

    T& cell(coord c)
    {
        WSS_ASSERT(c[0]>=0);
        WSS_ASSERT(c[0]<edge);
        WSS_ASSERT(c[1]>=0);
        WSS_ASSERT(c[1]<edge);

        return cells.get()[c[0]+c[1]*edge];
    }

private:
    int edge;
    std::unique_ptr<T[]> cells;
};

template<typename CellType, typename TextToCell>
std::optional<board<CellType>> make_board(
        std::vector<std::vector<char>> const& lines,
        TextToCell const& mapping)
{
    auto edge{ssize(lines)};
    board<CellType> result{int(edge)};

    for (auto row_index{0}; row_index!=edge; ++row_index) {
        auto const& line{lines[row_index]};
        auto const num_fields{ssize(line)};
        if (num_fields!=edge) {
            fmt::print(stderr,
                    "error: input row #{} has {} fields, expected {}.\n",
                    row_index+1, num_fields, edge);
            return std::nullopt;
        }

        for (auto column_index{0}; column_index!=edge; ++column_index) {
            auto const field{line[column_index]};
            auto const cell_found{mapping(field)};
            if (!cell_found) {
                fmt::print(
                        stderr,
                        "Unrecognised field, '{}', in row #{}, column #{}.\n",
						(char)field,
                        row_index+1, column_index+1);
                return std::nullopt;
            }

            result.cell(coord{column_index,row_index}) = *cell_found;
        }
    }

    return result;
}

template<typename CellType, typename TextToCell>
std::optional<board<CellType>> load_board(
        gsl::cstring_span<> filename,
        TextToCell const& mapping)
{
    auto const buffer{load_buffer(filename)};
    if (!buffer) {
        return std::nullopt;
    }

    auto const fields{parse_grid(*buffer)};

    return make_board<CellType>(fields, mapping);
}

template<typename T>
void transpose(board<T>& b)
{
    auto const edge = ssize(b);
    for (auto row = 0; row != edge; ++row)
    {
        for (auto column = 0; column != row; ++ column)
        {
            std::swap(b.cell(coord{column, row}), b.cell(coord{row, column}));
        }
    }
}

#endif //WSS_BOARD_H

max_stars_repo_name: johnmcfarlane/wss

max_issues_repo_name: johnmcfarlane/wss

max_forks_repo_name: johnmcfarlane/wss


text:
/* $Id$      */
/*--------------------------------------------------------------------*/
/*;  Copyright (C) 2008-2016                                          */
/*;  Associated Universities, Inc. Washington DC, USA.                */
/*;                                                                   */
/*;  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., 675 Massachusetts Ave, Cambridge,     */
/*;  MA 02139, USA.                                                   */
/*;                                                                   */
/*;Correspondence about this software should be addressed as follows: */
/*;         Internet email: bcotton@nrao.edu.                         */
/*;         Postal address: William Cotton                            */
/*;                         National Radio Astronomy Observatory      */
/*;                         520 Edgemont Road                         */
/*;                         Charlottesville, VA 22903-2475 USA        */
/*--------------------------------------------------------------------*/
#ifndef OBITSPECTRUMFIT_H
#define OBITSPECTRUMFIT_H

#include "Obit.h"
#include "ObitErr.h"
#include "ObitImage.h"
#include "ObitBeamShape.h"
#include "ObitThread.h"
#include "ObitInfoList.h"
#ifdef HAVE_GSL
#include <gsl/gsl_multifit_nlin.h>
#endif /* HAVE_GSL */

/*-------- Obit:  Merx mollis mortibus nuper ------------------*/
/**
 * \file ObitSpectrumFit.h
 *
 * ObitSpectrumFit Class for fitting spectra to image pixels
 *
 * This class does least squares fitting of log(s) as a polynomial in log($\nu$).
 * Either an image cube or a set of single plane images at arbitrary
 * frequencies may be fitted.
 * The result is an image cube of Log(S) with multiples of powers of log($\nu$)
 * as the planes.
 * The function ObitSpectrumFitEval will evaluate this fit and return an image
 * with the flux densities at the desired frequencies.
 *
 * \section ObitSpectrumFitaccess Creators and Destructors
 * An ObitSpectrumFit will usually be created using ObitSpectrumFitCreate which allows
 * specifying a name for the object as well as other information.
 *
 * A copy of a pointer to an ObitSpectrumFit should always be made using the
 * #ObitSpectrumFitRef function which updates the reference count in the object.
 * Then whenever freeing an ObitSpectrumFit or changing a pointer, the function
 * #ObitSpectrumFitUnref will decrement the reference count and destroy the object
 * when the reference count hits 0.
 * There is no explicit destructor.
 */

/*--------------Class definitions-------------------------------------*/
/** ObitSpectrumFit Class structure. */
typedef struct {
#include "ObitSpectrumFitDef.h"   /* this class definition */
} ObitSpectrumFit;

/*----------------- Macroes ---------------------------*/
/**
 * Macro to unreference (and possibly destroy) an ObitSpectrumFit
 * returns a ObitSpectrumFit*.
 * in = object to unreference
 */
#define ObitSpectrumFitUnref(in) ObitUnref (in)

/**
 * Macro to reference (update reference count) an ObitSpectrumFit.
 * returns a ObitSpectrumFit*.
 * in = object to reference
 */
#define ObitSpectrumFitRef(in) ObitRef (in)

/**
 * Macro to determine if an object is the member of this or a
 * derived class.
 * Returns TRUE if a member, else FALSE
 * in = object to reference
 */
#define ObitSpectrumFitIsA(in) ObitIsA (in, ObitSpectrumFitGetClass())

/*---------------Public functions---------------------------*/
/** Public: Class initializer. */
void ObitSpectrumFitClassInit (void);

/** Public: Default Constructor. */
ObitSpectrumFit* newObitSpectrumFit (gchar* name);

/** Public: Create/initialize ObitSpectrumFit structures */
ObitSpectrumFit* ObitSpectrumFitCreate (gchar* name, olong nterm);
/** Typedef for definition of class pointer structure */
typedef ObitSpectrumFit* (*ObitSpectrumFitCreateFP) (gchar* name,
						     olong nterm);

/** Public: ClassInfo pointer */
gconstpointer ObitSpectrumFitGetClass (void);

/** Public: Copy (deep) constructor. */
ObitSpectrumFit* ObitSpectrumFitCopy  (ObitSpectrumFit *in,
				       ObitSpectrumFit *out, ObitErr *err);

/** Public: Copy structure. */
void ObitSpectrumFitClone (ObitSpectrumFit *in, ObitSpectrumFit *out,
			   ObitErr *err);

/** Public: Fit spectrum to an image cube */
void ObitSpectrumFitCube (ObitSpectrumFit* in, ObitImage *inImage,
			  ObitImage *outImage, ObitErr *err);
/** Typedef for definition of class pointer structure */
typedef void(*ObitSpectrumFitCubeFP) (ObitSpectrumFit* in, ObitImage *inImage,
				      ObitImage *outImage, ObitErr *err);

/** Public: Fit spectrum to an array of images */
void ObitSpectrumFitImArr (ObitSpectrumFit* in, olong nimage, ObitImage **imArr,
			   ObitImage *outImage, ObitErr *err);
/** Typedef for definition of class pointer structure */
typedef void(*ObitSpectrumFitImArrFP) (ObitSpectrumFit* in, olong nimage, ObitImage **imArr,
			   ObitImage *outImage, ObitErr *err);

/* Do actual fitting */
void ObitSpectrumFitter (ObitSpectrumFit* in, ObitErr *err);
typedef void(*ObitSpectrumFitterFP) (ObitSpectrumFit* in, ObitErr *err);

/** Public: Evaluate spectrum */
void ObitSpectrumFitEval (ObitSpectrumFit* in, ObitImage *inImage,
			  odouble outFreq, ObitImage *outImage, ObitErr *err);
/** Typedef for definition of class pointer structure */
typedef void(*ObitSpectrumFitEvalFP) (ObitSpectrumFit* in, ObitImage *inImage,
				      odouble outFreq, ObitImage *outImage,
				      ObitErr *err);
/** Private: Write output image */
void ObitSpectrumWriteOutput (ObitSpectrumFit* in, ObitImage *outImage,
			      ObitErr *err);
typedef void (*ObitSpectrumWriteOutputFP) (ObitSpectrumFit* in, ObitImage *outImage,
					   ObitErr *err);
/** Public: Fit single spectrum */
ofloat* ObitSpectrumFitSingle (olong nfreq, olong nterm, odouble refFreq, odouble *freq,
			       ofloat *flux, ofloat *sigma, gboolean doBrokePow,
			       ObitErr *err);
/** Typedef for definition of class pointer structure */
typedef ofloat*(*ObitSpectrumFitSingleFP) (olong nfreq, olong nterm, odouble refFreq,
					   odouble *freq, ofloat *flux, ofloat *sigma,
					   gboolean doBrokePow, ObitErr *err);

/** Public: Make fitting arg structure */
gpointer ObitSpectrumFitMakeArg (olong nfreq, olong nterm,
				 odouble refFreq, odouble *freq,
				 gboolean doBrokePow,
				 ofloat **out, ObitErr *err);

/** Public: Fit single spectrum using arg */
void ObitSpectrumFitSingleArg (gpointer arg, ofloat *flux, ofloat  *sigma,
			       ofloat *out);

/** Public: Kill fitting arg structure */
void ObitSpectrumFitKillArg (gpointer arg);
/*----------- ClassInfo Structure -----------------------------------*/
/**
 * ClassInfo Structure.
 * Contains class name, a pointer to any parent class
 * (NULL if none) and function pointers.
 */
typedef struct  {
#include "ObitSpectrumFitClassDef.h"
} ObitSpectrumFitClassInfo;

#endif /* OBITFSPECTRUMFIT_H */

max_stars_repo_name: sarrvesh/Obit

max_issues_repo_name: sarrvesh/Obit

max_forks_repo_name: sarrvesh/Obit


text:
#include <stdio.h>
#include <string.h>
#include <gsl/gsl_statistics_double.h>
#include "tz_error.h"
#include "tz_constant.h"
#include "tz_image_lib.h"
#include "tz_objdetect.h"
#include "tz_imatrix.h"
#include "tz_stack_math.h"
#include "tz_stack_bwdist.h"
#include "tz_pixel_array.h"
#include "tz_stack_bwmorph.h"
#include "tz_stack_draw.h"
#include "tz_dmatrix.h"

INIT_EXCEPTION_MAIN(e)

int main(int argc, char* argv[])
{
  char neuron_name[100];
  if (argc == 1) {
    strcpy(neuron_name, "fly_neuron");
  } else {
    strcpy(neuron_name, argv[1]);
  }

  char file_path[100];

  sprintf(file_path, "../data/%s/objseed.tif", neuron_name);
  Stack *seedimg = Read_Stack(file_path);
  Object_3d_List *objs = Stack_Find_Object_N(seedimg, NULL, 1, 0, 26);

  Stack *stack2 = NULL;
  Object_3d *obj = NULL;
  stack2 = NULL;

  Object_3d_List *objs2 = NULL;

  sprintf(file_path, "../data/%s/bundle.tif", neuron_name);
  stack2 = Read_Stack(file_path);

  while (objs != NULL) {
    obj = Stack_Grow_Object(stack2, 1, objs->data->voxels[0]);
    if (obj != NULL) {
      if (obj->size > 1000) {
	if (objs2 == NULL) {
	  objs2 = Object_3d_List_New();
	  objs2->data = obj;
	} else {
	  Object_3d_List_Add(&objs2, obj);
	  Print_Object_3d_Info(obj);
	}
      } else {
	Kill_Object_3d(obj);
      }
    }
    objs = objs->next;
  }

  Print_Object_3d_List_Compact(objs2);

  Stack *bundle = Make_Stack(GREY, stack2->width, stack2->height, stack2->depth);
  Zero_Stack(bundle);

  Stack_Draw_Objects_Bw(bundle, objs2, 1);

  sprintf(file_path, "../data/%s/grow_bundle.tif", neuron_name);
  Write_Stack(file_path, bundle);

  return 0;
}


text:
//
//  geru.h
//  Linear Algebra Template Library
//
//  Created by Rodney James on 12/23/11.
//  Copyright (c) 2011 University of Colorado Denver. All rights reserved.
//

#ifndef _geru_h
#define _geru_h

/// @file geru.h Performs complex vector outer product.

#include "latl.h"

namespace LATL
{
   /// @brief Performs a vector outer product of two complex vectors.
   ///
   /// For a complex matrix A, complex vectors x and y, and complex scalar alpha,
   ///
   ///        A := alpha * x * y.' + A
   ///
   /// is computed.
   /// @return 0 if success.
   /// @return -i if the ith argument is invalid.
   /// @tparam real_t Floating point type.
   /// @param m Specifies the number of rows of the matrix A.  m>=0
   /// @param n Specifies the number of coumns of the matrix A.  n>=0
   /// @param alpha Complex scalar.
   /// @param x Pointer to complex vector x.
   /// @param incx Increment of the vector x.  x!=0
   /// @param y Pointer to complex vector y.
   /// @param incy Increment of the vector y.  y!=0
   /// @param A Pointer to complex m-by-n matrix A.
   /// @param ldA Column length of matrix A.  ldA>=m.
   /// @ingroup BLAS

   template <typename real_t>
   int GERU(int_t m, int_t n, complex<real_t> alpha, complex<real_t> *x, int_t incx, complex<real_t> *y, int_t incy, complex<real_t> *A, int_t ldA)
   {
      using std::conj;

      const complex<real_t> zero(0.0,0.0);
      int_t i,j,kx,jy,ix;

      if(m<0)
         return -1;
      else if(n<0)
         return -2;
      else if(incx==0)
         return -5;
      else if(incy==0)
         return -7;
      else if(ldA<m)
         return -9;
      else if((m==0)||(n==0)||(alpha==zero))
         return 0;

      jy=(incy>0)?0:(1-n)*incy;
      kx=(incx>0)?0:(1-m)*incx;

      if(incx==1)
      {
         for(j=0;j<n;j++)
         {
            for(i=0;i<m;i++)
               A[i]+=x[i]*alpha*y[jy];
            jy+=incy;
            A+=ldA;
         }
      }
      else
      {
         for(j=0;j<n;j++)
         {
            ix=kx;
            for(i=0;i<m;i++)
            {
               A[i]+=x[ix]*alpha*y[jy];
               ix+=incx;
            }
            A+=ldA;
            jy+=incy;
         }
      }
      return 0;
   }
#ifdef __latl_cblas
#include <cblas.h>

   template <> int GERU<float>(int_t m, int_t n, complex<float> alpha, complex<float> *x, int_t incx, complex<float> *y, int_t incy, complex<float> *A, int_t ldA)
   {
      if(m<0)
         return -1;
      else if(n<0)
         return -2;
      else if(incx==0)
         return -5;
      else if(incy==0)
         return -7;
      else if(ldA<m)
         return -9;

      cblas_cgeru(CblasColMajor,m,n,&alpha,x,incx,y,incy,A,ldA);

      return 0;
   }

   template <> int GERU<double>(int_t m, int_t n, complex<double> alpha, complex<double> *x, int_t incx, complex<double> *y, int_t incy, complex<double> *A, int_t ldA)
   {
      if(m<0)
         return -1;
      else if(n<0)
         return -2;
      else if(incx==0)
         return -5;
      else if(incy==0)
         return -7;
      else if(ldA<m)
         return -9;

      cblas_cgeru(CblasColMajor,m,n,&alpha,x,incx,y,incy,A,ldA);

      return 0;
   }

#endif
}
#endif

max_stars_repo_name: langou/latl

max_issues_repo_name: langou/latl

max_forks_repo_name: langou/latl


text:
#pragma once
#include "GradUtil.h"

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

#include <limits>
#include <math.h>
#include <vector>
#include "BasicError.h"
#include "Util.h"

using namespace std;

class ValueGrad;
class DistanceGrad {
public:
	double dist;
	gsl_vector* grad;
	bool set; // TODO: this bit is not necessary if the default dist is 0
	DistanceGrad(double d, gsl_vector* g): dist(d), grad(g), set(true) {}
	~DistanceGrad(void) {
		gsl_vector_free(grad);
	}
	string print() {
		stringstream str;
		if (set) {
			str << "Dist: " << dist;
		} else {
			str << "Dist: NOT SET";
		}
		return str.str();
	}
	string printFull() {
		stringstream str;
		if (set) {
			str << "Dist: " << dist << endl;
			str << "DGrads: ";
			for (int i = 0; i < grad->size; i++) {
				str << gsl_vector_get(grad, i) << ", ";
			}
            if (dist > 1e5 || gsl_blas_dnrm2(grad) > 1e5) {
                str << "LARGE VALUES" << endl;
            }
		} else {
			str << "Dist: NOT SET";
		}
		return str.str();
	}

	static void dg_and(DistanceGrad* m, DistanceGrad* f, DistanceGrad* d);
	static void dg_or(DistanceGrad* m, DistanceGrad* f, DistanceGrad* d);
	static void dg_not(DistanceGrad* m, DistanceGrad* d);
	static void dg_ite(DistanceGrad* b, DistanceGrad* m, DistanceGrad* f, DistanceGrad* d);

	static double dg_copy(DistanceGrad* m, gsl_vector* grad);
	static void dg_copy(DistanceGrad* m, DistanceGrad* o);

	static double dg_times(DistanceGrad* m, DistanceGrad* f, gsl_vector* grad);

    static void dg_ite(DistanceGrad* m, DistanceGrad* f, double dval, gsl_vector* dgrad, DistanceGrad* o);

    static bool same(DistanceGrad* m, DistanceGrad* f);

    static double dg_combine(DistanceGrad* m, DistanceGrad* f);
    static void dg_combine(DistanceGrad* m, DistanceGrad* f, DistanceGrad* o);

    static double dg_combine(DistanceGrad* m, DistanceGrad* f, double cval, gsl_vector* grad, int bv, gsl_vector* o);
    static double dg_combine(DistanceGrad* m, DistanceGrad* f, double cval, int bv);

    static void dg_combine(vector<DistanceGrad*>& cdists, DistanceGrad* vdist, vector<ValueGrad*>& cvals, int bv, DistanceGrad* o);
};

max_stars_repo_name: natebragg/sketch-backend

max_issues_repo_name: natebragg/sketch-backend

max_forks_repo_name: natebragg/sketch-backend


text:
/*
 * Author : Pierre Schnizer
 * Date   : February 2005
 *
 */
#include <pygsl/error_helpers.h>
#include <pygsl/block_helpers.h>
#include <pygsl/pygsl_features.h>
#include <string.h>

#include <gsl/gsl_fft.h>
#include <gsl/gsl_fft_complex.h>
#include <gsl/gsl_fft_real.h>
#include <gsl/gsl_fft_halfcomplex.h>
#include <gsl/gsl_fft_complex_float.h>
#include <gsl/gsl_fft_real_float.h>
#include <gsl/gsl_fft_halfcomplex_float.h>
#include <gsl/gsl_blas.h>
#ifdef _PYGSL_GSL_HAS_WAVELET
#define forward forward_wavelet
#define backward backward_wavelet
#include <gsl/gsl_wavelet.h>
#include <gsl/gsl_wavelet2d.h>
#undef forward
#undef backward
#else
/*
#error Need GSL version 1.6 or higher with wavelet implementation! If you need it for a lower version
#error add the different conditional compilation!
*/
#endif /* gsl version > 1.5*/
static PyObject *module = NULL;
static const char filename[] = __FILE__;

/* The types */
#include "transformtypes.h"
/* The documentation strings */
#include "doc.h"
/* The python type which allows to allocate and pass the work space */
#include "space.c"
/* Function to copy the different arrays. */
#include "arraycopy.c"
/*
 * Macros and functions which set the information up and call the
 * PyGSL_transform functions
 */
#include "fft.c"
/*
 * The wavelet type and the methods to call the different transform functions.
 */
#ifdef _PYGSL_GSL_HAS_WAVELET
#include "wavelet.c"
#define PyGSL_WAVELET_TRANSFORM_TYPE(name) \
       {#name,             PyGSL_wavelet_init_ ## name,              METH_VARARGS, NULL}, \
       {#name "_centered", PyGSL_wavelet_init_ ## name ## _centered, METH_VARARGS, NULL},
#else
#define PyGSL_WAVELET_TRANSFORM_TYPE(name)
#endif
/*
 * The real workers.
 */
#include "core.c"

static PyMethodDef transformMethods[] = {
	/* workspace init functions  */
	{"complex_workspace",           PyGSL_transform_space_init_COMPLEX_WORKSPACE,            METH_VARARGS, (char*)cws_doc},
	{"complex_wavetable",           PyGSL_transform_space_init_COMPLEX_WAVETABLE,            METH_VARARGS, (char*)cwt_doc},
	{"real_workspace",              PyGSL_transform_space_init_REAL_WORKSPACE,               METH_VARARGS, (char*)rws_doc},
	{"real_wavetable",              PyGSL_transform_space_init_REAL_WAVETABLE,               METH_VARARGS, (char*)rwt_doc},
	{"halfcomplex_wavetable",       PyGSL_transform_space_init_HALFCOMPLEX_WAVETABLE,        METH_VARARGS, (char*)hcwt_doc},
	{"complex_workspace_float",     PyGSL_transform_space_init_COMPLEX_WORKSPACE_FLOAT,      METH_VARARGS, (char*)cws_doc},
	{"complex_wavetable_float",     PyGSL_transform_space_init_COMPLEX_WAVETABLE_FLOAT,      METH_VARARGS, (char*)cwt_doc},
	{"real_workspace_float",        PyGSL_transform_space_init_REAL_WORKSPACE_FLOAT,         METH_VARARGS, (char*)rws_doc},
	{"real_wavetable_float",        PyGSL_transform_space_init_REAL_WAVETABLE_FLOAT,         METH_VARARGS, (char*)rwt_doc},
	{"halfcomplex_wavetable_float", PyGSL_transform_space_init_HALFCOMPLEX_WAVETABLE_FLOAT,  METH_VARARGS, (char*)hcwt_doc},
#ifdef _PYGSL_GSL_HAS_WAVELET
	{"wavelet_workspace",           PyGSL_transform_space_init_WAVELET_WORKSPACE,            METH_VARARGS, (char*)ww_doc},
#endif
	/* transform functions */
	PyGSL_TRANSFORM_FD_FUNCTION("complex_forward",             fft_complex_forward,              cf_doc)
	PyGSL_TRANSFORM_FD_FUNCTION("complex_backward",            fft_complex_backward,             cb_doc)
	PyGSL_TRANSFORM_FD_FUNCTION("complex_inverse",             fft_complex_inverse,              ci_doc)
	PyGSL_TRANSFORM_FD_FUNCTION("complex_radix2_forward",      fft_complex_radix2_forward,       cf_doc_r2)
	PyGSL_TRANSFORM_FD_FUNCTION("complex_radix2_backward",     fft_complex_radix2_backward,      cb_doc_r2)
	PyGSL_TRANSFORM_FD_FUNCTION("complex_radix2_inverse",      fft_complex_radix2_inverse,       ci_doc_r2)
	PyGSL_TRANSFORM_FD_FUNCTION("complex_radix2_dif_forward",  fft_complex_radix2_dif_forward,   cf_doc_r2_dif)
	PyGSL_TRANSFORM_FD_FUNCTION("complex_radix2_dif_backward", fft_complex_radix2_dif_backward,  cb_doc_r2_dif)
	PyGSL_TRANSFORM_FD_FUNCTION("complex_radix2_dif_inverse",  fft_complex_radix2_dif_inverse,   ci_doc_r2_dif)
	PyGSL_TRANSFORM_FD_FUNCTION("real_transform",              fft_real_transform,               rt_doc)
	PyGSL_TRANSFORM_FD_FUNCTION("halfcomplex_transform",       fft_halfcomplex_transform,        hc_doc)
	PyGSL_TRANSFORM_FD_FUNCTION("halfcomplex_inverse",         fft_halfcomplex_inverse,          hi_doc)
	PyGSL_TRANSFORM_FD_FUNCTION("real_radix2_transform",       fft_real_radix2_transform,        rt_doc_r2)
	PyGSL_TRANSFORM_FD_FUNCTION("halfcomplex_radix2_transform",fft_halfcomplex_radix2_transform, hc_doc_r2)
	PyGSL_TRANSFORM_FD_FUNCTION("halfcomplex_radix2_inverse",  fft_halfcomplex_radix2_inverse,   hi_doc_r2)
	/* helper functions */
	{"halfcomplex_radix2_unpack",         PyGSL_fft_halfcomplex_radix2_unpack,       METH_VARARGS, (char*)un_doc_r2},
	{"halfcomplex_radix2_unpack_float",   PyGSL_fft_halfcomplex_radix2_unpack_float, METH_VARARGS, (char*)float_doc},
	/* wavelet inits */
#ifdef _PYGSL_GSL_HAS_WAVELET
	PyGSL_WAVELET_TRANSFORM_TYPE(daubechies)
	PyGSL_WAVELET_TRANSFORM_TYPE(haar)
	PyGSL_WAVELET_TRANSFORM_TYPE(bspline)
#endif
	{NULL, NULL, 0, NULL} /* Sentinel */
};

/*
 * Set the various function pointers for the different transforms. See the
 * structure _pygsl_transform_func_s for the functions. Some architectures do
 * not allow to initalise function pointers in static structures on the heap.
 * (some solaris versions? )
 */
#define PYGSL_INIT_FUNCS(helpers, space, table, spacet, tablet) \
	helpers.space_alloc = (pygsl_transform_helpn_t *) gsl_fft_ ## space ## _alloc; \
	helpers.space_free  = (pygsl_transform_help_t *)  gsl_fft_ ## space ## _free;  \
	helpers.table_alloc = (pygsl_transform_helpn_t *) gsl_fft_ ## table ## _alloc; \
	helpers.table_free  = (pygsl_transform_help_t *)  gsl_fft_ ## table ## _free; \
	helpers.space_type  = spacet; \
	helpers.table_type  = tablet; \

#define PYGSL_INIT_FUNCS_DF(helpers, space, table, spacet, tablet) \
       PYGSL_INIT_FUNCS(helpers ## _funcs, space, table, spacet, tablet) \
       PYGSL_INIT_FUNCS(helpers ## _float ## _funcs, space ## _float, table ## _float, spacet ## _FLOAT, tablet ## _FLOAT)

static void
init_helpers(void)
{
	FUNC_MESS_BEGIN();
	PYGSL_INIT_FUNCS_DF(complex,     complex_workspace, complex_wavetable, COMPLEX_WORKSPACE, COMPLEX_WAVETABLE)
	PYGSL_INIT_FUNCS_DF(real,        real_workspace,    real_wavetable, REAL_WORKSPACE, REAL_WAVETABLE)
	PYGSL_INIT_FUNCS_DF(halfcomplex, real_workspace,    halfcomplex_wavetable, REAL_WORKSPACE, HALFCOMPLEX_WAVETABLE)

        DEBUG_MESS(3, "PyArray_FLOAT   = %d ", PyArray_FLOAT  );
	DEBUG_MESS(3, "PyArray_DOUBLE  = %d ", PyArray_DOUBLE );
	DEBUG_MESS(3, "PyArray_CFLOAT  = %d ", PyArray_CFLOAT );
	DEBUG_MESS(3, "PyArray_CDOUBLE = %d ", PyArray_CDOUBLE);

#ifdef _PYGSL_GSL_HAS_WAVELET
	DEBUG_MESS(4, "%s @ %p", "daubechies", gsl_wavelet_daubechies);
	DEBUG_MESS(4, "%s @ %p", "daubechies_centered", gsl_wavelet_daubechies_centered);
	DEBUG_MESS(4, "%s @ %p", "haar", gsl_wavelet_haar);
	DEBUG_MESS(4, "%s @ %p", "haar_centered", gsl_wavelet_haar_centered);
	DEBUG_MESS(4, "%s @ %p", "bspline", gsl_wavelet_bspline);
	DEBUG_MESS(4, "%s @ %p", "bspline_centered", gsl_wavelet_bspline_centered);
#endif
	FUNC_MESS_END();

}

DL_EXPORT(void) init_transform(void)
{
     	PyObject *m = NULL, *dict = NULL, *item = NULL;

	FUNC_MESS_BEGIN();
	PyGSL_transform_space_pytype.ob_type = &PyType_Type;
#ifdef _PYGSL_GSL_HAS_WAVELET
	PyGSL_wavelet_pytype.ob_type = &PyType_Type;
#endif
	m = Py_InitModule("_transform", transformMethods);
	module = m;
	init_pygsl();
	init_helpers();
	if (m == NULL)
		return;

	dict = PyModule_GetDict(m);
	if (dict == NULL)
		return;

	if (!(item = PyString_FromString(transform_module_doc))){
		PyErr_SetString(PyExc_ImportError,
				"I could not generate module doc string!");
		return;
	}
	if (PyDict_SetItemString(dict, "__doc__", item) != 0){
		PyErr_SetString(PyExc_ImportError,
				"I could not init doc string!");
		return;
	}
	FUNC_MESS_END();
	return;
}

/*
 * Local Variables:
 * mode: C
 * c-file-style: "python"
 * End:
 */

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


text:
//
// Created by Harold on 2020/9/30.
//

#ifndef M_MATH_M_INTEGRATION_H
#define M_MATH_M_INTEGRATION_H

#include <gsl/gsl_integration.h>

namespace M_MATH {
    class IntegrationQAG {
    public:
        typedef double (*func_type)(double x, void* params);
    public:
        IntegrationQAG() : IntegrationQAG(1000, 1e-6) { }

        IntegrationQAG(size_t workspace_sz, double accuracy)
                : workspace_size(workspace_sz), accuracy_bound(accuracy), F() {
            workspace = gsl_integration_workspace_alloc(workspace_size);
        }

        ~IntegrationQAG() {
            gsl_integration_workspace_free(workspace);
        }

        void bind_func(func_type f, void* params) {
            F.function = f;
            F.params = params;
        }

        void integrate(double x0, double x1, double* result, double* error) {
            // default use key = 1: GSL_INTEG_GAUSS15
            gsl_integration_qag(&F, x0, x1, 0, accuracy_bound, workspace_size, 1, workspace, result, error);
        }
    private:
        gsl_function F;
        size_t workspace_size;
        double accuracy_bound;
        gsl_integration_workspace* workspace;
    };

    // change lambda function to function ptr
    struct Lambda {
        template<typename Tret, typename T>
        static Tret lambda_ptr_exec(double x, void* data) {
            return (Tret) (*(T*)fn<T>())(x, data);
        }

        template<typename Tret = void, typename Tfp = Tret(*)(double, void*), typename T>
        static Tfp ptr(T& t) {
            fn<T>(&t);
            return (Tfp) lambda_ptr_exec<Tret, T>;
        }

        template<typename T>
        static void* fn(void* new_fn = nullptr) {
            static void* fn;
            if (new_fn != nullptr)
                fn = new_fn;
            return fn;
        }
    };
}

#endif //M_MATH_M_INTEGRATION_H


text:
/* multirobust.c
 *
 * Copyright (C) 2013 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.
 *
 * This module contains routines related to robust linear least squares. The
 * algorithm used closely follows the publications:
 *
 * [1] DuMouchel, W. and F. O'Brien (1989), "Integrating a robust
 * option into a multiple regression computing environment,"
 * Computer Science and Statistics:  Proceedings of the 21st
 * Symposium on the Interface, American Statistical Association
 *
 * [2] Street, J.O., R.J. Carroll, and D. Ruppert (1988), "A note on
 * computing robust regression estimates via iteratively
 * reweighted least squares," The American Statistician, v. 42,
 * pp. 152-154.
 */

//#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_sort_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>

#include "linear_common.c"

static int robust_test_convergence(const gsl_vector *c_prev, const gsl_vector *c,
                                   const double tol);
static double robust_madsigma(const gsl_vector *x, gsl_multifit_robust_workspace *w);
static double robust_robsigma(const gsl_vector *r, const double s,
                              const double tune, gsl_multifit_robust_workspace *w);
static double robust_sigma(const double s_ols, const double s_rob,
                           gsl_multifit_robust_workspace *w);
static int robust_covariance(const double sigma, gsl_matrix *cov,
                             gsl_multifit_robust_workspace *w);

/*
gsl_multifit_robust_alloc
  Allocate a robust workspace

Inputs: T - robust weighting algorithm
        n - number of observations
        p - number of model parameters

Return: pointer to workspace
*/

gsl_multifit_robust_workspace *
gsl_multifit_robust_alloc(const gsl_multifit_robust_type *T,
                          const size_t n, const size_t p)
{
  gsl_multifit_robust_workspace *w;

  if (n < p)
    {
      GSL_ERROR_VAL("observations n must be >= p", GSL_EINVAL, 0);
    }

  w = calloc(1, sizeof(gsl_multifit_robust_workspace));
  if (w == 0)
    {
      GSL_ERROR_VAL("failed to allocate space for multifit_robust struct",
                    GSL_ENOMEM, 0);
    }

  w->n = n;
  w->p = p;
  w->type = T;
  /* bdavis */
  //w->maxiter = 100; /* maximum iterations */
  w->maxiter = 5; /* maximum iterations */
  /* bdavis */
  w->tune = w->type->tuning_default;

  w->multifit_p = gsl_multifit_linear_alloc(n, p);
  if (w->multifit_p == 0)
    {
      GSL_ERROR_VAL("failed to allocate space for multifit_linear struct",
                    GSL_ENOMEM, 0);
    }

  w->r = gsl_vector_alloc(n);
  if (w->r == 0)
    {
      GSL_ERROR_VAL("failed to allocate space for residuals",
                    GSL_ENOMEM, 0);
    }

  w->weights = gsl_vector_alloc(n);
  if (w->weights == 0)
    {
      GSL_ERROR_VAL("failed to allocate space for weights", GSL_ENOMEM, 0);
    }

  w->c_prev = gsl_vector_alloc(p);
  if (w->c_prev == 0)
    {
      GSL_ERROR_VAL("failed to allocate space for c_prev", GSL_ENOMEM, 0);
    }

  w->resfac = gsl_vector_alloc(n);
  if (w->resfac == 0)
    {
      GSL_ERROR_VAL("failed to allocate space for residual factors",
                    GSL_ENOMEM, 0);
    }

  w->psi = gsl_vector_alloc(n);
  if (w->psi == 0)
    {
      GSL_ERROR_VAL("failed to allocate space for psi", GSL_ENOMEM, 0);
    }

  w->dpsi = gsl_vector_alloc(n);
  if (w->dpsi == 0)
    {
      GSL_ERROR_VAL("failed to allocate space for dpsi", GSL_ENOMEM, 0);
    }

  w->QSI = gsl_matrix_alloc(p, p);
  if (w->QSI == 0)
    {
      GSL_ERROR_VAL("failed to allocate space for QSI", GSL_ENOMEM, 0);
    }

  w->D = gsl_vector_alloc(p);
  if (w->D == 0)
    {
      GSL_ERROR_VAL("failed to allocate space for D", GSL_ENOMEM, 0);
    }

  w->workn = gsl_vector_alloc(n);
  if (w->workn == 0)
    {
      GSL_ERROR_VAL("failed to allocate space for workn", GSL_ENOMEM, 0);
    }

  w->stats.sigma_ols = 0.0;
  w->stats.sigma_mad = 0.0;
  w->stats.sigma_rob = 0.0;
  w->stats.sigma = 0.0;
  w->stats.Rsq = 0.0;
  w->stats.adj_Rsq = 0.0;
  w->stats.rmse = 0.0;
  w->stats.sse = 0.0;
  w->stats.dof = n - p;
  w->stats.weights = w->weights;
  w->stats.r = w->r;

  return w;
} /* gsl_multifit_robust_alloc() */

/*
gsl_multifit_robust_free()
  Free memory associated with robust workspace
*/

void
gsl_multifit_robust_free(gsl_multifit_robust_workspace *w)
{
  if (w->multifit_p)
    gsl_multifit_linear_free(w->multifit_p);

  if (w->r)
    gsl_vector_free(w->r);

  if (w->weights)
    gsl_vector_free(w->weights);

  if (w->c_prev)
    gsl_vector_free(w->c_prev);

  if (w->resfac)
    gsl_vector_free(w->resfac);

  if (w->psi)
    gsl_vector_free(w->psi);

  if (w->dpsi)
    gsl_vector_free(w->dpsi);

  if (w->QSI)
    gsl_matrix_free(w->QSI);

  if (w->D)
    gsl_vector_free(w->D);

  if (w->workn)
    gsl_vector_free(w->workn);

  free(w);
} /* gsl_multifit_robust_free() */

int
gsl_multifit_robust_tune(const double tune, gsl_multifit_robust_workspace *w)
{
  w->tune = tune;
  return GSL_SUCCESS;
}

const char *
gsl_multifit_robust_name(const gsl_multifit_robust_workspace *w)
{
  return w->type->name;
}

gsl_multifit_robust_stats
gsl_multifit_robust_statistics(const gsl_multifit_robust_workspace *w)
{
  return w->stats;
}

/*
gsl_multifit_robust()
  Perform robust iteratively reweighted linear least squares
fit

Inputs: X     - design matrix of basis functions
        y     - right hand side vector
        c     - (output) model coefficients
        cov   - (output) covariance matrix
        w     - workspace
*/

int
gsl_multifit_robust(const gsl_matrix * X,
                    const gsl_vector * y,
                    gsl_vector * c,
                    gsl_matrix * cov,
                    gsl_multifit_robust_workspace *w)
{
  /* check matrix and vector sizes */
  if (X->size1 != y->size)
    {
      GSL_ERROR
        ("number of observations in y does not match rows of matrix X",
         GSL_EBADLEN);
    }
  else if (X->size2 != c->size)
    {
      GSL_ERROR ("number of parameters c does not match columns of matrix X",
                 GSL_EBADLEN);
    }
  else if (cov->size1 != cov->size2)
    {
      GSL_ERROR ("covariance matrix is not square", GSL_ENOTSQR);
    }
  else if (c->size != cov->size1)
    {
      GSL_ERROR
        ("number of parameters does not match size of covariance matrix",
         GSL_EBADLEN);
    }
  else if (X->size1 != w->n || X->size2 != w->p)
    {
      GSL_ERROR
        ("size of workspace does not match size of observation matrix",
         GSL_EBADLEN);
    }
  else
    {
      int s;
      double chisq;
      const double tol = GSL_SQRT_DBL_EPSILON;
      int converged = 0;
      size_t numit = 0;
      const size_t n = y->size;
      double sigy = gsl_stats_sd(y->data, y->stride, n);
      double sig_lower;
      size_t i;

      /*
       * if the initial fit is very good, then finding outliers by comparing
       * them to the residual standard deviation is difficult. Therefore we
       * set a lower bound on the standard deviation estimate that is a small
       * fraction of the standard deviation of the data values
       */
      sig_lower = 1.0e-6 * sigy;
      if (sig_lower == 0.0)
        sig_lower = 1.0;

      /* compute initial estimates using ordinary least squares */
      s = gsl_multifit_linear(X, y, c, cov, &chisq, w->multifit_p);
      if (s)
        return s;

      /* save Q S^{-1} of original matrix */
      gsl_matrix_memcpy(w->QSI, w->multifit_p->QSI);
      gsl_vector_memcpy(w->D, w->multifit_p->D);

      /* compute statistical leverage of each data point */
      s = gsl_linalg_SV_leverage(w->multifit_p->A, w->resfac);
      if (s)
        return s;

      /* correct residuals with factor 1 / sqrt(1 - h) */
      for (i = 0; i < n; ++i)
        {
          double h = gsl_vector_get(w->resfac, i);

          if (h > 0.9999)
            h = 0.9999;

          gsl_vector_set(w->resfac, i, 1.0 / sqrt(1.0 - h));
        }

      /* compute residuals from OLS fit r = y - X c */
      s = gsl_multifit_linear_residuals(X, y, c, w->r);
      if (s)
        return s;

      /* compute estimate of sigma from ordinary least squares */
      w->stats.sigma_ols = gsl_blas_dnrm2(w->r) / sqrt((double) w->stats.dof);

      while (!converged && ++numit <= w->maxiter)
        {
          double sig;

          /* adjust residuals by statistical leverage (see DuMouchel and O'Brien) */
          s = gsl_vector_mul(w->r, w->resfac);
          if (s)
            return s;

          /* compute estimate of standard deviation using MAD */
          sig = robust_madsigma(w->r, w);

          /* scale residuals by standard deviation and tuning parameter */
          gsl_vector_scale(w->r, 1.0 / (GSL_MAX(sig, sig_lower) * w->tune));

          /* compute weights using these residuals */
          s = w->type->wfun(w->r, w->weights);
          if (s)
            return s;

          gsl_vector_memcpy(w->c_prev, c);

          /* solve weighted least squares with new weights */
          s = gsl_multifit_wlinear(X, w->weights, y, c, cov, &chisq, w->multifit_p);
          if (s)
            return s;

          /* compute new residuals r = y - X c */
          s = gsl_multifit_linear_residuals(X, y, c, w->r);
          if (s)
            return s;

          converged = robust_test_convergence(w->c_prev, c, tol);
        }

      /* compute final MAD sigma */
      w->stats.sigma_mad = robust_madsigma(w->r, w);

      /* compute robust estimate of sigma */
      w->stats.sigma_rob = robust_robsigma(w->r, w->stats.sigma_mad, w->tune, w);

      /* compute final estimate of sigma */
      w->stats.sigma = robust_sigma(w->stats.sigma_ols, w->stats.sigma_rob, w);

      /* store number of iterations */
      w->stats.numit = numit;

      {
        double dof = (double) w->stats.dof;
        double rnorm = w->stats.sigma * sqrt(dof); /* see DuMouchel, sec 4.2 */
        double ss_err = rnorm * rnorm;
        double ss_tot = gsl_stats_tss(y->data, y->stride, n);

        /* compute R^2 */
        w->stats.Rsq = 1.0 - ss_err / ss_tot;

        /* compute adjusted R^2 */
        w->stats.adj_Rsq = 1.0 - (1.0 - w->stats.Rsq) * (n - 1.0) / dof;

        /* compute rmse */
        w->stats.rmse = sqrt(ss_err / dof);

        /* store SSE */
        w->stats.sse = ss_err;
      }

      /* calculate covariance matrix = sigma^2 (X^T X)^{-1} */
      s = robust_covariance(w->stats.sigma, cov, w);
      if (s)
        return s;

      /* raise an error if not converged */
/* bdavis */
      /* Eliminating this check is to avoid an error when iterations */
      /* exceed 5.  A better solution is probably recommended, such as */
      /* reverting to default of 100 if an input specification is not */
      /* enabled and defined. */
      /* bdavis@usgs.gov */
      /*
      if (numit > w->maxiter)
        {
          GSL_ERROR("maximum iterations exceeded", GSL_EMAXITER);
        }
      */
/* bdavis */

      return s;
    }
} /* gsl_multifit_robust() */

/* Estimation of values for given x */
int
gsl_multifit_robust_est(const gsl_vector * x, const gsl_vector * c,
                        const gsl_matrix * cov, double *y, double *y_err)
{
  int s = gsl_multifit_linear_est(x, c, cov, y, y_err);

  return s;
}

/***********************************
 * INTERNAL ROUTINES               *
 ***********************************/

/*
robust_test_convergence()
  Test for convergence in robust least squares

Convergence criteria:

|c_i^(k) - c_i^(k-1)| <= tol * max(|c_i^(k)|, |c_i^(k-1)|)

for all i. k refers to iteration number.

Inputs: c_prev - coefficients from previous iteration
        c      - coefficients from current iteration
        tol    - tolerance

Return: 1 if converged, 0 if not
*/

static int
robust_test_convergence(const gsl_vector *c_prev, const gsl_vector *c,
                        const double tol)
{
  size_t p = c->size;
  size_t i;

  for (i = 0; i < p; ++i)
    {
      double ai = gsl_vector_get(c_prev, i);
      double bi = gsl_vector_get(c, i);

      if (fabs(bi - ai) > tol * GSL_MAX(fabs(ai), fabs(bi)))
        return 0; /* not yet converged */
    }

  /* converged */
  return 1;
} /* robust_test_convergence() */

/*
robust_madsigma()
  Estimate the standard deviation of the residuals using
the Median-Absolute-Deviation (MAD) of the residuals,
throwing away the smallest p residuals.

See: Street et al, 1988

Inputs: r - vector of residuals
        w - workspace
*/

static double
robust_madsigma(const gsl_vector *r, gsl_multifit_robust_workspace *w)
{
  gsl_vector_view v;
  double sigma;
  size_t n = r->size;
  const size_t p = w->p;
  size_t i;

  /* copy |r| into workn */
  for (i = 0; i < n; ++i)
    {
      gsl_vector_set(w->workn, i, fabs(gsl_vector_get(r, i)));
    }

  gsl_sort_vector(w->workn);

  /*
   * ignore the smallest p residuals when computing the median
   * (see Street et al 1988)
   */
  v = gsl_vector_subvector(w->workn, p - 1, n - p + 1);
  sigma = gsl_stats_median_from_sorted_data(v.vector.data, v.vector.stride, v.vector.size) / 0.6745;

  return sigma;
} /* robust_madsigma() */

/*
robust_robsigma()
  Compute robust estimate of sigma so that
sigma^2 * inv(X' * X) is a reasonable estimate of
the covariance for robust regression. Based heavily
on the equations of Street et al, 1988.

Inputs: r    - vector of residuals y - X c
        s    - sigma estimate using MAD
        tune - tuning constant
        w    - workspace
*/

static double
robust_robsigma(const gsl_vector *r, const double s,
                const double tune, gsl_multifit_robust_workspace *w)
{
  double sigma;
  size_t i;
  const size_t n = w->n;
  const size_t p = w->p;
  const double st = s * tune;
  double a, b, lambda;

  /* compute u = r / sqrt(1 - h) / st */
  gsl_vector_memcpy(w->workn, r);
  gsl_vector_mul(w->workn, w->resfac);
  gsl_vector_scale(w->workn, 1.0 / st);

  /* compute w(u) and psi'(u) */
  w->type->wfun(w->workn, w->psi);
  w->type->psi_deriv(w->workn, w->dpsi);

  /* compute psi(u) = u*w(u) */
  gsl_vector_mul(w->psi, w->workn);

  /* Street et al, Eq (3) */
  a = gsl_stats_mean(w->dpsi->data, w->dpsi->stride, n);

  /* Street et al, Eq (5) */
  b = 0.0;
  for (i = 0; i < n; ++i)
    {
      double psi_i = gsl_vector_get(w->psi, i);
      double resfac = gsl_vector_get(w->resfac, i);
      double fac = 1.0 / (resfac*resfac); /* 1 - h */

      b += fac * psi_i * psi_i;
    }
  b /= (double) (n - p);

  /* Street et al, Eq (5) */
  lambda = 1.0 + ((double)p)/((double)n) * (1.0 - a) / a;

  sigma = lambda * sqrt(b) * st / a;

  return sigma;
} /* robust_robsigma() */

/*
robust_sigma()
  Compute final estimate of residual standard deviation, using
the OLS and robust sigma estimates.

This equation is taken from DuMouchel and O'Brien, sec 4.1:
\hat{\sigma_R}

Inputs: s_ols - OLS sigma
        s_rob - robust sigma
        w     - workspace

Return: final estimate of sigma
*/

static double
robust_sigma(const double s_ols, const double s_rob,
             gsl_multifit_robust_workspace *w)
{
  double sigma;
  const size_t p = w->p;
  const size_t n = w->n;

  /* see DuMouchel and O'Brien, sec 4.1 */
  sigma = GSL_MAX(s_rob,
                  sqrt((s_ols*s_ols*p*p + s_rob*s_rob*n) /
                       (p*p + n)));

  return sigma;
} /* robust_sigma() */

/*
robust_covariance()
  Calculate final covariance matrix, defined as:

  sigma * (X^T X)^{-1}

Inputs: sigma - residual standard deviation
        cov   - (output) covariance matrix
        w     - workspace
*/

static int
robust_covariance(const double sigma, gsl_matrix *cov,
                  gsl_multifit_robust_workspace *w)
{
  int s = 0;
  const size_t p = w->p;
  const double s2 = sigma * sigma;
  size_t i, j;
  gsl_matrix *QSI = w->QSI;
  gsl_vector *D = w->D;

  /* Form variance-covariance matrix cov = s2 * (Q S^-1) (Q S^-1)^T */

  for (i = 0; i < p; i++)
    {
      gsl_vector_view row_i = gsl_matrix_row (QSI, i);
      double d_i = gsl_vector_get (D, i);

      for (j = i; j < p; j++)
        {
          gsl_vector_view row_j = gsl_matrix_row (QSI, j);
          double d_j = gsl_vector_get (D, j);
          double s;

          gsl_blas_ddot (&row_i.vector, &row_j.vector, &s);

          gsl_matrix_set (cov, i, j, s * s2 / (d_i * d_j));
          gsl_matrix_set (cov, j, i, s * s2 / (d_i * d_j));
        }
    }

  return s;
} /* robust_covariance() */

max_stars_repo_name: agroimpacts/imager

max_issues_repo_name: agroimpacts/imager

max_forks_repo_name: agroimpacts/imager