text:
#pragma once
#include "parser.h"
#include "format.h"
#include "nameutils.h"
#include "utils.h"
#include <sfun/string_utils.h>
#include <cmdlime/errors.h>
#include <gsl/gsl>
#include <algorithm>
#include <sstream>
#include <iomanip>
#include <functional>
#include <optional>

namespace cmdlime::detail{
namespace str = sfun::string_utils;

template <FormatType formatType>
class GNUParser : public Parser<formatType>
{
    using Parser<formatType>::Parser;
    using FindMode = typename Parser<formatType>::FindMode;

    void preProcess() override
    {
        checkNames();
        foundParam_.clear();
        foundParamPrefix_.clear();
    }

    void process(const std::string& token) override
    {
        if (str::startsWith(token, "--") && token.size() > 2)
            processCommand(token);
        else if (str::startsWith(token, "-") && token.size() > 1)
            processShortCommand(token);
        else if (!foundParam_.empty()){
            this->readParam(foundParam_, token);
            foundParam_.clear();
        }
        else
            this->readArg(token);
    }

    void postProcess() override
    {
        if (!foundParam_.empty())
            throw ParsingError{"Parameter '" + foundParamPrefix_ + foundParam_ + "' value can't be empty"};
    }

    void processCommand(std::string command)
    {
        command = str::after(command, "--");
        auto paramValue = std::optional<std::string>{};
        if (command.find('=') != std::string::npos){
            paramValue = str::after(command, "=");
            command = str::before(command, "=");
        }
        if (isParamOrFlag(command) &&
            !foundParam_.empty() &&
            this->readMode_ != Parser<formatType>::ReadMode::ExitFlagsAndCommands)
            throw ParsingError{"Parameter '" + foundParamPrefix_ + foundParam_ + "' value can't be empty"};
        if (this->findParam(command, FindMode::Name) || this->findParamList(command, FindMode::Name)){
            if (paramValue.has_value())
                this->readParam(command, paramValue.value());
            else{
                foundParam_ = command;
                foundParamPrefix_ = "--";
            }
        }
        else if (this->findFlag(command, FindMode::Name))
            this->readFlag(command);
        else if (this->readMode_ != Parser<formatType>::ReadMode::ExitFlagsAndCommands)
            throw ParsingError{"Encountered unknown parameter or flag '--" + command + "'"};
    }

    void processShortCommand(std::string command)
    {
        auto possibleNumberArg = command;
        command = str::after(command, "-");
        if (isShortParamOrFlag(command)){
            if (!foundParam_.empty() && this->readMode_ != Parser<formatType>::ReadMode::ExitFlagsAndCommands)
                throw ParsingError{"Parameter '" + foundParamPrefix_ + foundParam_ + "' value can't be empty"};
            parseShortCommand(command);
        }
        else if (isNumber(possibleNumberArg))
            this->readArg(possibleNumberArg);
        else if (this->readMode_ != Parser<formatType>::ReadMode::ExitFlagsAndCommands)
            throw ParsingError{"Encountered unknown parameter or flag '-" + command + "'"};
    }

    void parseShortCommand(const std::string& command)
    {
        if (command.empty())
            throw ParsingError{"Flags and parameters must have a name"};
        auto paramValue = std::string{};
        for(auto ch : command){
            auto opt = std::string{ch};
            if (!foundParam_.empty())
                paramValue += opt;
            else if (this->findFlag(opt, FindMode::ShortName))
                this->readFlag(opt);
            else if (this->findParam(opt, FindMode::ShortName)){
                foundParam_ = opt;
                foundParamPrefix_ = "-";
            }
            else if (this->findParamList(opt, FindMode::ShortName)){
                foundParam_ = opt;
                foundParamPrefix_ = "-";
            }
            else if (this->readMode_ != Parser<formatType>::ReadMode::ExitFlagsAndCommands)
                throw ParsingError{"Unknown option '" + opt + "' in command '-" + command + "'"};
        }
        if (!foundParam_.empty() && !paramValue.empty()){
            this->readParam(foundParam_, paramValue);
            foundParam_.clear();
        }
    }

    void checkLongNames()
    {
        auto check = [](const OptionInfo& var, const std::string& varType){
            if (!std::isalpha(var.name().front()))
                throw ConfigError{varType + "'s name '" + var.name() + "' must start with an alphabet character"};
            if (var.name().size() > 1){
                auto nonSupportedCharIt = std::find_if(var.name().begin() + 1, var.name().end(), [](char ch){return !std::isalnum(ch) && ch != '-';});
                if (nonSupportedCharIt != var.name().end())
                    throw ConfigError{varType + "'s name '" + var.name() + "' must consist of alphanumeric characters and hyphens"};
            }
        };
        this->forEachParamInfo([check](const OptionInfo& var){
            check(var, "Parameter");
        });
        this->forEachParamListInfo([check](const OptionInfo& var){
            check(var, "Parameter");
        });
        this->forEachFlagInfo([check](const OptionInfo& var){
            check(var, "Flag");
        });
    }

    void checkShortNames()
    {
        auto check = [](const OptionInfo& var, const std::string& varType){
            if (var.shortName().empty())
                return;
            if (var.shortName().size() != 1)
                throw ConfigError{varType + "'s short name '" + var.shortName() + "' can't have more than one symbol"};
            if (!std::isalnum(var.shortName().front()))
                throw ConfigError{varType + "'s short name '" + var.shortName() + "' must be an alphanumeric character"};
        };
        this->forEachParamInfo([check](const OptionInfo& var){
            check(var, "Parameter");
        });
        this->forEachParamListInfo([check](const OptionInfo& var){
            check(var, "Parameter");
        });
        this->forEachFlagInfo([check](const OptionInfo& var){
            check(var, "Flag");
        });
    }

    void checkNames()
    {
        checkLongNames();
        checkShortNames();
    }

    bool isParamOrFlag(const std::string& str)
    {
        if (str.empty())
            return false;
        return this->findFlag(str, FindMode::Name) ||
               this->findParam(str, FindMode::Name) ||
               this->findParamList(str, FindMode::Name);
    }

    bool isShortParamOrFlag(const std::string& str)
    {
        if (str.empty())
            return false;
        auto opt = str.substr(0,1);
        return this->findFlag(opt, FindMode::ShortName) ||
               this->findParam(opt, FindMode::ShortName) ||
               this->findParamList(opt, FindMode::ShortName);
    }

private:
    std::string foundParam_;
    std::string foundParamPrefix_;
};

class GNUNameProvider{
public:
    static std::string name(const std::string& optionName)
    {
        Expects(!optionName.empty());
        return toKebabCase(optionName);
    }

    static std::string shortName(const std::string& optionName)
    {
        Expects(!optionName.empty());
        return toLowerCase(optionName.substr(0, 1));
    }

    static std::string fullName(const std::string& optionName)
    {
        Expects(!optionName.empty());
        return toKebabCase(optionName);
    }

    static std::string valueName(const std::string& typeName)
    {
        Expects(!typeName.empty());
        return toKebabCase(templateType(typeNameWithoutNamespace(typeName)));
    }
};

class GNUOutputFormatter{
public:
    static std::string paramUsageName(const IParam& param)
    {
        auto stream = std::stringstream{};
        if (param.isOptional())
            stream << "[" << paramPrefix() << param.info().name() << " <" << param.info().valueName() << ">]";
        else
            stream << paramPrefix() << param.info().name() << " <" << param.info().valueName() << ">";
        return stream.str();
    }

    static std::string paramListUsageName(const IParamList& param)
    {
        auto stream = std::stringstream{};
        if (param.isOptional())
            stream << "[" << paramPrefix() << param.info().name() << " <" << param.info().valueName() << ">...]";
        else
            stream << paramPrefix() << param.info().name() << " <" << param.info().valueName() << ">...";
        return stream.str();
    }

    static std::string paramDescriptionName(const IParam& param, int indent = 0)
    {
        auto stream = std::stringstream{};
        stream << std::setw(indent);
        if (!param.info().shortName().empty())
            stream << "-" << param.info().shortName() << ", ";
        else
            stream << " " << "   ";
        stream << "--" << param.info().name() << " <" << param.info().valueName() << ">";
        return stream.str();
    }

    static std::string paramListDescriptionName(const IParamList& param, int indent = 0)
    {
        auto stream = std::stringstream{};
        stream << std::setw(indent);
        if (!param.info().shortName().empty())
            stream << "-" << param.info().shortName() << ", ";
        else
            stream << " " << "   ";
        stream << "--" << param.info().name() << " <" << param.info().valueName() << ">";
        return stream.str();
    }

    static std::string paramPrefix()
    {
        return "--";
    }

    static std::string flagUsageName(const IFlag& flag)
    {
        auto stream = std::stringstream{};
        stream << "[" << flagPrefix() << flag.info().name() << "]";
        return stream.str();
    }

    static std::string flagDescriptionName(const IFlag& flag, int indent = 0)
    {
        auto stream = std::stringstream{};
        stream << std::setw(indent) ;
        if (!flag.info().shortName().empty())
            stream << "-" << flag.info().shortName() << ", ";
        else
            stream << " " << "   ";
        stream << "--" << flag.info().name();
        return stream.str();
    }

    static std::string flagPrefix()
    {
        return "--";
    }

    static std::string argUsageName(const IArg& arg)
    {
        auto stream = std::stringstream{};
        stream << "<" << arg.info().name() << ">";
        return stream.str();
    }

    static std::string argDescriptionName(const IArg& arg, int indent = 0)
    {
        auto stream = std::stringstream{};
        if (indent)
            stream << std::setw(indent) << " ";
        stream << "<" << arg.info().name() << "> (" << arg.info().valueName() << ")";
        return stream.str();
    }

    static std::string argListUsageName(const IArgList& argList)
    {
        auto stream = std::stringstream{};
        if (argList.isOptional())
            stream << "[" << argList.info().name() << "...]";
        else
            stream << "<" << argList.info().name() << "...>";
        return stream.str();
    }

    static std::string argListDescriptionName(const IArgList& argList, int indent = 0)
    {
        auto stream = std::stringstream{};
        if (indent)
            stream << std::setw(indent) << " ";
        stream << "<" << argList.info().name() << "> (" << argList.info().valueName() << ")";
        return stream.str();
    }

};

template<>
struct Format<FormatType::GNU>
{
    using parser = GNUParser<FormatType::GNU>;
    using nameProvider = GNUNameProvider;
    using outputFormatter = GNUOutputFormatter;
    static constexpr bool shortNamesEnabled = true;
};

}

ext: h

max_stars_repo_name: GerHobbelt/hypertextcpp

max_issues_repo_name: GerHobbelt/hypertextcpp

max_forks_repo_name: GerHobbelt/hypertextcpp


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

#include "allvars.h"
#include "proto.h"

void init()
{
  int i;
  double dT, T1, T2;

  strcpy(fname_results, "data/results.txt");
  fp_results = fopen(fname_results, "w");

#ifdef JAVELIN
  fprintf(fp_results, "JAVELIN-single top-hat transfer function.\n");
  printf("JAVELIN.\n");
  nc_lim_low = nc_lim_up = 1;
#elif defined TOPHAT
  fprintf(fp_results, "multiple top-hat transfer function.\n");
  printf("TOPHATs.\n");
#else
  fprintf(fp_results, "multiple Gaussian transfer function.\n");
  printf("GAUSSIANs.\n");
#endif

  memory_alloc_data();
  read_data();

  nall = ncon + nline;
  nall_data = ncon_data + nline_data;
  memory_alloc();

  scale_con = scale_line = 1.0;
  scale_light_curves();

  memcpy(Fall_data, Fcon_data, ncon_data*sizeof(double));
  memcpy(Fall_data+ncon_data, Fline_data, nline_data*sizeof(double));

  len_con = Tcon_data[ncon_data-1] - Tcon_data[0];
  len_line = Tline_data[nline_data-1] - Tline_data[0];

  cad_con = 0.0;
  for(i=0; i<ncon_data-1; i++)
  {
    cad_con += Tcon_data[i+1] - Tcon_data[i];
  }
  cad_con /= (ncon_data-1);

  cad_line = 0.0;
  for(i=0; i<nline_data-1; i++)
  {
    cad_line += Tline_data[i+1] - Tline_data[i];
  }
  cad_line /= (nline_data-1);

  printf("len: %f %f\n", len_con, len_line);
  printf("cad: %f %f\n", cad_con, cad_line);

  T1 = Tcon_data[0];
  T2 = Tcon_data[ncon_data-1];

  dT = T2-T1;
  T1 -= 0.1*dT;
  T2 += 0.1*dT;
  dT = (T2-T1)/(ncon-1.0);
  for(i=0; i<ncon; i++)
    Tcon[i] = T1 + dT *i;

  T1 = Tline_data[0];
  T2 = Tline_data[nline_data-1];

  dT = T2-T1;
  T1 -= 0.1*dT;
  T2 += 0.1*dT;
  dT = (T2-T1)/(nline-1.0);
  for(i=0; i<nline; i++)
    Tline[i] = T1 + dT *i;
}

void scale_light_curves()
{
  int i;
  double mean, norm;

  mean = 0.0;
  norm = 0.0;
  for(i=0; i<ncon_data; i++)
  {
    mean += Fcon_data[i]/(Fcerrs_data[i] * Fcerrs_data[i]);
    norm += 1.0/(Fcerrs_data[i] * Fcerrs_data[i]);
  }
  mean /= norm;
  scale_con = mean;

  for(i=0; i<ncon_data;i++)
  {
    Fcon_data[i] /= mean;
    Fcerrs_data[i] /= mean;
  }

  mean = 0.0;
  norm = 0.0;
  for(i=0; i<nline_data; i++)
  {
    mean += Fline_data[i]/(Flerrs_data[i]*Flerrs_data[i]);
    norm += 1.0/(Flerrs_data[i]*Flerrs_data[i]);
  }
  mean /= norm;
  scale_line = mean;
  for(i=0; i<nline_data;i++)
  {
    Fline_data[i] /= mean;
    Flerrs_data[i] /= mean;
  }

  printf("scale: %e %e\n", scale_con, scale_line);
  fprintf(fp_results, "scale: %e %e\n", scale_con, scale_line);
}

/*
 *
 */
void memory_alloc_data()
{
  Tcon_data = malloc(ndata_max*sizeof(double));
  if(Tcon_data==NULL)
  {
    strcpy(str_error_exit, "Tcon_data");
    error_exit(7);
  }

  Fcon_data = malloc(ndata_max*sizeof(double));
  if(Fcon_data==NULL)
  {
    strcpy(str_error_exit, "Fcon_data");
    error_exit(7);
  }

  Fcerrs_data = malloc(ndata_max*sizeof(double));
  if(Fcerrs_data==NULL)
  {
    strcpy(str_error_exit, "Fcerrs_data");
    error_exit(7);
  }

  Tline_data = malloc(ndata_max*sizeof(double));
  if(Tline_data==NULL)
  {
    strcpy(str_error_exit, "Tline_data");
    error_exit(7);
  }

  Fline_data = malloc(ndata_max*sizeof(double));
  if(Fline_data==NULL)
  {
    strcpy(str_error_exit, "Fline_data");
    error_exit(7);
  }

  Flerrs_data = malloc(ndata_max*sizeof(double));
  if(Flerrs_data==NULL)
  {
    strcpy(str_error_exit, "Flerrs_data");
    error_exit(7);
  }
}
void memory_alloc()
{
  workspace_ipiv = malloc(ndata_max*sizeof(double));
  workspace = malloc(10*ndata_max*sizeof(double) + 10);

  Fall_data = array_malloc(nall_data);

  Smat = array_malloc(nall_data*nall_data);
  Nmat = array_malloc(nall_data*nall_data);
  Cmat = array_malloc(nall_data*nall_data);
  ICmat = array_malloc(nall_data*nall_data);
  ICvmat = array_malloc(nall_data*nall_data);

  USmat = array_malloc(nall*nall_data);
  ASmat = array_malloc(nall*nall);

  Tmat1 = array_malloc(nall_data*nall_data);
  Tmat2 = array_malloc(nall_data*nall_data);
  Tmat3 = array_malloc(nall_data*nall_data);
  Tmat4 = array_malloc(nall_data*nall_data);

  Tcon = array_malloc(ncon);
  Fcon = array_malloc(ncon);
  Fcerrs = array_malloc(ncon);
  Tline = array_malloc(nline);
  Fline = array_malloc(nline);
  Flerrs = array_malloc(nline);

  theta_range = matrix_malloc(ntheta_max, 2);
  theta_best = array_malloc(ntheta_max);
  theta_best_var = array_malloc(ntheta_max*2);
  theta_input = array_malloc(ntheta_max);
  sigma_input = array_malloc(ntheta_max);
  theta_fixed = malloc(ntheta_max*sizeof(int));

  theta_best_con = array_malloc(ntheta_max);
  theta_best_var_con = array_malloc(ntheta_max*2);

  cov_matrix = array_malloc(ntheta_max * ntheta_max);

  grid_tau = array_malloc(nc_lim_up);
  TF_tau = array_malloc(ntau);
  TF = array_malloc(ntau);

  gsl_T = gsl_rng_default;
  gsl_r = gsl_rng_alloc (gsl_T);
  gsl_rng_set(gsl_r, time(NULL));
}

ext: c

max_stars_repo_name: LiyrAstroph/MICA

max_issues_repo_name: LiyrAstroph/MICA

max_forks_repo_name: LiyrAstroph/MICA


text:
#ifndef TTT_UTIL_TEXTURECACHE_H
#define TTT_UTIL_TEXTURECACHE_H

#include <SDL_render.h>

#include <gsl/pointers>
#include <gsl/util>

#include <vector>
#include <memory>
#include <string_view>

namespace ttt::gfx
{
	// A class for caching textures for reuse.
	class TextureCache final
	{
	public:
		using SizeType = gsl::index;

		SDL_Texture *load(std::string_view path, gsl::not_null<SDL_Renderer *> renderer);
		void unload(SizeType index) noexcept;
		void erase(SizeType index) noexcept;
		void clear() noexcept;

		[[nodiscard]] SDL_Texture *at(SizeType index) const noexcept;
		[[nodiscard]] SDL_Texture *operator [](SizeType index) const noexcept;

		[[nodiscard]] SizeType getSize() const noexcept;
		[[nodiscard]] SizeType getMaxSize() const noexcept;

	private:
		std::vector<std::unique_ptr<SDL_Texture, decltype(&SDL_DestroyTexture)>> textures;
	};

	inline void TextureCache::unload(SizeType index) noexcept { gsl::at(textures, index).reset(); }
	inline void TextureCache::erase(SizeType index) noexcept { textures.erase(textures.cbegin() + index); }
	inline void TextureCache::clear() noexcept { textures.clear(); }

	inline SDL_Texture *TextureCache::at(SizeType index) const noexcept { return gsl::at(textures, index).get(); }
	inline SDL_Texture *TextureCache::operator [](SizeType index) const noexcept { return gsl::at(textures, index).get(); }

	inline TextureCache::SizeType TextureCache::getSize() const noexcept { return textures.size(); }
	inline TextureCache::SizeType TextureCache::getMaxSize() const noexcept { return textures.max_size(); }
}

#endif

max_stars_repo_name: itsArtem/TicTacToe

max_issues_repo_name: itsArtem/TicTacToe

max_forks_repo_name: itsArtem/TicTacToe


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

#pragma once

#ifndef vdir_b62819a0_88a2_4eb7_9a92_6d287a152304_h
#define vdir_b62819a0_88a2_4eb7_9a92_6d287a152304_h

#include <gslib/config.h>
#include <gslib/tree.h>
#include <gslib/pool.h>

__gslib_begin__

class dirop
{
public:
    virtual void get_current_dir(string& str);
    virtual bool change_dir(const gchar* path);
    virtual bool make_dir(const gchar* path);
    virtual bool remove_dir(const gchar* path);
};

struct vdirnode;
struct vdirfile;
typedef vdirnode vdirfolder;
typedef _treenode_wrapper<vdirnode> vdirwrapper;
typedef tree<vdirnode, vdirwrapper> vdirtree;
typedef vdirtree::iterator vdiriter;
typedef vdirtree::const_iterator vdirciter;

struct vdirnode
{
    enum
    {
        dt_folder,
        dt_file,
        dt_tag,
    };
    string          name;

public:
    vdirnode() {}
    vdirnode(const gchar* n) { name.assign(n); }
    virtual uint get_tag() const { return dt_folder; }
    virtual const gchar* get_name() const { return name.c_str(); }
    void set_name(const gchar* str, int len) { name.assign(str, len); }
    void set_name(const string& str) { name = str; }
    void tracing() const;
};

struct vdirfile:
    public vdirnode
{
    vessel          vsl;

public:
    vdirfile() {}
    vdirfile(const gchar* n): vdirnode(n) {}
    virtual uint get_tag() const override { return dt_file; }

public:
    const gchar* get_postfix() const;
    void set_file_data(const void* ptr, int size);
    void set_file_data(const gchar* filename);
};

struct vdirtag:
    public vdirnode
{
    void*           binding;

public:
    vdirtag() { binding = 0; }
    vdirtag(const gchar* n): vdirnode(n) { binding = 0; }
    virtual uint get_tag() const override { return dt_tag; }
    void set_binding(void* p) { binding = p; }
    void* get_binding() const { return binding; }
};

class vdir:
    public vdirtree
{
public:
    vdiriter create_folder(vdiriter pos, const gchar* name);
    vdiriter create_file(vdiriter pos, const gchar* name);
    vdiriter create_tag(vdiriter pos, const gchar* name);
    void save(const gchar* path);
};

class vdirop:
    public dirop
{
public:
    vdirop(vdir& d): _dir(d) {}
    virtual void get_current_dir(string& str) override;
    virtual bool change_dir(const gchar* path) override;
    virtual bool make_dir(const gchar* path) override;
    virtual bool remove_dir(const gchar* path) override;

protected:
    vdir&           _dir;
    vdiriter        _curr;

public:
    vdiriter step_or_create(vdiriter pos, const gchar* name);
    vdiriter step_once(vdiriter pos, const gchar* name);
    vdiriter step_to(vdiriter pos, const gchar* path);
    vdiriter locate_folder(vdiriter pos, const gchar* name);
    vdiriter locate_file(vdiriter pos, const gchar* name);
    vdiriter create_folder(const gchar* name) { return _dir.create_folder(_curr, name); }
    vdiriter create_file(const gchar* name) { return _dir.create_file(_curr, name); }
    vdiriter get_current_iter() const { return _curr; }
    void rewind_curr() { _curr = _dir.get_root(); }
};

__gslib_end__

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
#include <fftw.h>
#include <stdio.h>
#include <stdlib.h>
#define N 16
main()
{
     int i;
     fftw_complex in[N], out[N];
     float fact;
     fftw_plan p,p2;
printf(" stuff in data\n");
     for (i=0;i<N;i++) {
     	in[i].re=(i*i);
     	in[i].im=1;
     }
printf(" create plans\n");
     p = fftw_create_plan(N, FFTW_FORWARD,  FFTW_ESTIMATE);
     p2= fftw_create_plan(N, FFTW_BACKWARD, FFTW_ESTIMATE);
printf(" do it\n");
     fftw_one(p, in, out);
     for (i=0;i<N;i++) {
	printf("%12.4f %12.4f\n",out[i].re,out[i].im);
     }
printf(" \n");
printf(" undo it\n");
     fftw_one(p2, out,in);
     fact=1.0/(float)N;
     for (i=0;i<N;i++) {
	printf("%10.2f %10.2f\n",in[i].re*fact,in[i].im*fact);
     }
printf(" clean up\n");
     fftw_destroy_plan(p);
     fftw_destroy_plan(p2);
}

max_stars_repo_name: timkphd/examples

max_issues_repo_name: timkphd/examples

max_forks_repo_name: timkphd/examples


text:
// lasp_math_raw.c
//
// last-edit-by: J.A. de Jong
//
// Description:
// Operations working on raw arrays of floating point numbers
//////////////////////////////////////////////////////////////////////
#define TRACERPLUS (-5)
#include "lasp_math_raw.h"
#if LASP_USE_BLAS
#include <cblas.h>
#endif

void d_elem_prod_d(d res[],
                   const d arr1[],
                   const d arr2[],
                   const us size) {

    #if LASP_USE_BLAS == 1

    #if LASP_DEBUG

    if(arr1 == arr2) {
        DBGWARN("d_elem_prod_d: Array 1 and array 2 point to the same"
                " memory. This results in pointer aliasing, for which"
                " testing is still to be done. Results might be"
                " unrealiable.");
    }

    #endif

    #if LASP_DOUBLE_PRECISION
    #define elem_prod_fun cblas_dsbmv
    #else
    #define elem_prod_fun cblas_ssbmv
    #endif
    /* These parameters do not matter for this specific case */
    const CBLAS_ORDER  mat_order= CblasColMajor;
    const CBLAS_UPLO   uplo = CblasLower;

    /* Extra multiplication factor */
    const d alpha = 1.0;

    /* void cblas_dsbmv(OPENBLAS_CONST enum CBLAS_ORDER order, */
    /*                  OPENBLAS_CONST enum CBLAS_UPLO Uplo, */
    /*                  OPENBLAS_CONST blasint N, */
    /*                  OPENBLAS_CONST blasint K, */
    /*                  OPENBLAS_CONST double alpha, */
    /*                  OPENBLAS_CONST double *A, */
    /*                  OPENBLAS_CONST blasint lda, */
    /*                  OPENBLAS_CONST double *X, */
    /*                  OPENBLAS_CONST blasint incX, */
    /*                  OPENBLAS_CONST double beta, */
    /*                  double *Y, */
    /*                  OPENBLAS_CONST blasint incY); */

    elem_prod_fun(mat_order,
                  uplo,
                  (blasint) size,
                  0,             // Just the diagonal; 0 super-diagonal bands
                  alpha,        /* Multiplication factor alpha */
                  arr1,
                  1,            /* LDA */
                  arr2,         /* x */
                  1, /* incX = 1 */
                  0.0,          /* Beta */
                  res,    /* The Y matrix to write to */
                  1); /* incY */
    #undef elem_prod_fun

    #else  /* No blas routines, routine is very simple, but here we
            * go! */
    DBGWARN("Performing slow non-blas vector-vector multiplication");
    for(us i=0;i<size;i++) {
        res[i] = arr1[i]*arr2[i];
    }
    #endif
}

void c_hadamard(c res[],
                const c arr1[],
                const c arr2[],
                const us size) {

    fsTRACE(15);
    uVARTRACE(15,size);
    dbgassert(arr1 && arr2 && res,NULLPTRDEREF);

    #if LASP_USE_BLAS == 1

    #if LASP_DEBUG

    if(arr1 == arr2) {
        DBGWARN("c_elem_prod_c: Array 1 and array 2 point to the same"
                " memory. This results in pointer aliasing, for which"
                " testing is still to be done. Results might be"
                " unrealiable.");
    }

    #endif  /* LASP_DEBUG */

    #if LASP_DOUBLE_PRECISION
    #define elem_prod_fun cblas_zgbmv
    #else
    #define elem_prod_fun cblas_cgbmv
    #endif

    c alpha = 1.0;
    c beta = 0.0;

    TRACE(15,"Calling " annestr(elem_prod_fun));
    uVARTRACE(15,size);
    elem_prod_fun(CblasColMajor,
                  CblasNoTrans,
                  (blasint) size, /* M: Number of rows */
                  (blasint) size, /* B: Number of columns */
                  0,              /* KL: Number of sub-diagonals */
                  0,              /* KU: Number of super-diagonals */
                  (d*) &alpha,    /* Multiplication factor */
                  (d*) arr2,      /* A */
                  1,              /* LDA */
                  (d*) arr1,      /* x */
                  1,              /* incX = 1 */
                  (d*) &beta,
                  (d*) res,    /* The Y matrix to write to */
                  1);          /* incY (increment in res) */

    #undef elem_prod_fun

    #else  /* No blas routines, routine is very simple, but here we
            * go! */
    DBGWARN("Performing slower non-blas vector-vector multiplication");
    for(us i=0;i<size;i++) {
        res[i] = arr1[i]*arr2[i];
    }
    #endif
    feTRACE(15);
}

//////////////////////////////////////////////////////////////////////

max_stars_repo_name: asceenl/lasp

max_issues_repo_name: asceenl/lasp

max_forks_repo_name: asceenl/lasp


text:
/**
 * @file modalbeamformer.h
 * @brief Beamforming in the spherical harmonics domain.
 * @author Kenichi Kumatani
 */
#ifndef MODALBEAMFORMER_H
#define MODALBEAMFORMER_H

#include <stdio.h>
#include <assert.h>
#include <float.h>

#include <gsl/gsl_block.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_fft_complex.h>
#include <common/refcount.h>
#include "common/jexception.h"

#include "stream/stream.h"
//#include "stream/pyStream.h"
#include "beamformer/spectralinfoarray.h"
#include "modulated/modulated.h"
#include "beamformer/beamformer.h"

// ----- definition for class `ModeAmplitudeCalculator' -----
//
gsl_complex modeAmplitude(int order, double ka);

class ModeAmplitudeCalculator {
 public:
  ModeAmplitudeCalculator(int order, float minKa=0.01, float maxKa=20, float wid=0.01);

  ~ModeAmplitudeCalculator();

  gsl_vector_complex *get() const { return mode_amplitude_; }

private:
  gsl_vector_complex *mode_amplitude_;
  float minKa_;
  float maxKa_;
  float wid_;
};

typedef refcount_ptr<ModeAmplitudeCalculator> ModeAmplitudeCalculatorPtr;

// ----- definition for class `EigenBeamformer' -----
//

/**
   @class EigenBeamformer
   @brief This beamformer is implemented based on Meyer and Elko's ICASSP paper.
          In Boaz Rafaely's paper, this method is referred to as the phase-mode beamformer
   @usage
   1) construct this object, bf = EigenBeamformer(...)
   2) set the radious of the spherical array with bf.set_eigenmike_geometry() or bf.set_array_geometry(..).
   3) set the look direction with bf.set_look_direction(..)
   4) process each block with bf.next() until it hits the end
*/
class EigenBeamformer : public  SubbandDS {
public:
  EigenBeamformer( unsigned sampleRate, unsigned fftLen = 512, bool halfBandShift = false, unsigned NC=1, unsigned maxOrder=4, bool normalizeWeight=false, const String& nm = "EigenBeamformer");
  ~EigenBeamformer();
  virtual const gsl_vector_complex* next(int frame_no = -5);
  virtual void reset();
  virtual unsigned dim() const { return dim_;}

  void set_sigma2(float sigma2){ sigma2_ = sigma2; }
  void set_weight_gain(float wgain){ wgain_ = wgain; }
  void set_eigenmike_geometry();
  void set_array_geometry(double a,  gsl_vector *theta_s, gsl_vector *phi_s);
  virtual void set_look_direction(double theta, double phi);

  const gsl_matrix_complex *mode_amplitudes();
  const gsl_vector *array_geometry(int type); // type==0 -> theta, type==1 -> phi
  virtual gsl_matrix *beampattern(unsigned fbinX, double theta = 0, double phi = 0,
                                  double minTheta=-M_PI, double maxTheta=M_PI,
                                  double minPhi=-M_PI, double maxPhi=M_PI,
                                  double widthTheta=0.1, double widthPhi=0.1 );
  /**
     @brief obtain the spherical transformation coefficients at each frame
     @return spherical harmonics transformation coefficients at the current frame
   */
  virtual SnapShotArrayPtr snapshot_array() const { return(st_snapshot_array_); }
  virtual SnapShotArrayPtr snapshot_array2() const { return(snapshot_array_); }

  const gsl_matrix_complex *blocking_matrix(unsigned fbinX, unsigned unitX=0 ) const {
    return (bfweight_vec_[unitX]->B())[fbinX];
  }

#ifdef ENABLE_LEGACY_BTK_API
  void setSigma2(float sigma2){ set_sigma2(sigma2); }
  void setWeightGain(float wgain){ set_weight_gain(wgain); }
  void setEigenMikeGeometry(){ set_eigenmike_geometry(); }
  void setArrayGeometry(double a,  gsl_vector *theta_s, gsl_vector *phi_s){ set_array_geometry(a, theta_s, phi_s); }
  virtual void setLookDirection(double theta, double phi){ set_look_direction(theta, phi); }
  const gsl_matrix_complex *getModeAmplitudes(){ return mode_amplitudes(); }
  const gsl_vector *getArrayGeometry(int type){ return array_geometry(type);}
  virtual gsl_matrix *getBeamPattern( unsigned fbinX, double theta = 0, double phi = 0,
                                      double minTheta=-M_PI, double maxTheta=M_PI,
                                      double minPhi=-M_PI, double maxPhi=M_PI,
                                      double widthTheta=0.1, double widthPhi=0.1 ){
    return beampattern(fbinX, theta, phi, minTheta, maxTheta, minPhi, maxPhi, widthTheta, widthPhi);
  }
  virtual SnapShotArrayPtr getSnapShotArray(){ return snapshot_array(); }
  virtual SnapShotArrayPtr getSnapShotArray2(){ return snapshot_array2(); }
  const gsl_matrix_complex *getBlockingMatrix(unsigned fbinX, unsigned unitX=0){
    return blocking_matrix(fbinX, unitX);
  }
#endif

 protected:
  virtual void calc_weights_( unsigned fbinX, gsl_vector_complex *weights );
  virtual bool calc_spherical_harmonics_at_each_position_( gsl_vector *theta_s, gsl_vector *phi_s ); // need to be tested!!
  virtual bool calc_steering_unit_(  int unitX=0, bool isGSC=false );
  virtual bool alloc_steering_unit_( int unitN=1 );
  void alloc_image_( bool flag=true );
  bool calc_mode_amplitudes_();

  unsigned samplerate_;
  unsigned NC_;
  unsigned maxOrder_;
  unsigned dim_; // the number of the spherical harmonics transformation coefficients
  bool weights_normalized_;
  gsl_matrix_complex *mode_mplitudes_; // [maxOrder_] the mode amplitudes.
  gsl_vector_complex *F_; // Spherical Transform coefficients [dim_]
  gsl_vector_complex **sh_s_; // Conjugate of spherical harmonics at each sensor position [dim_][nChan]: Y_n^{m*}
  SnapShotArrayPtr st_snapshot_array_; // for compatibility with a post-filtering object

  double theta_; // look direction
  double phi_;   // look direction

  double a_;     // the radius of the rigid sphere.
  gsl_vector *theta_s_; // sensor positions
  gsl_vector *phi_s_;   // sensor positions
  gsl_matrix *beampattern_;
  gsl_vector *WNG_; // white noise gain
  float wgain_; //
  float sigma2_; // dialog loading
};

typedef Inherit<EigenBeamformer,  SubbandDSPtr> EigenBeamformerPtr;

// ----- definition for class DOAEstimatorSRPEB' -----
//

/**
   @class DOAEstimatorSRPEB
   @brief estimate the direction of arrival based on the maximum steered response power
   @usage
   1) construct this object, doaEstimator = DOAEstimatorSRPEB(...)
   2) set the radious of the spherical array, doaEstimator.set_eigenmike_geometry() or doaEstimator.set_array_geometry(..).
   3) process each block, doaEstimator.next()
   4) get the N-best hypotheses at the current instantaneous frame through doaEstimator.nbest_doas()
   5) do doaEstimator.getFinalNBestHypotheses() after a static segment is processed.
      You can then obtain the averaged N-best hypotheses of the static segment with doaEstimator.nbest_doas().
*/
class DOAEstimatorSRPEB :
  public DOAEstimatorSRPBase, public EigenBeamformer {
  public:
  DOAEstimatorSRPEB( unsigned nBest, unsigned sampleRate, unsigned fftLen = 512, bool halfBandShift = false, unsigned NC=1, unsigned maxOrder=4, bool normalizeWeight=false, const String& nm = "DirectionEstimatorSRPBase");
  ~DOAEstimatorSRPEB();

  const gsl_vector_complex* next(int frame_no = -5);
  void reset();

protected:
  virtual void  calc_steering_unit_table_();
  virtual float calc_response_power_( unsigned uttX );
};

typedef Inherit<DOAEstimatorSRPEB, EigenBeamformerPtr> DOAEstimatorSRPEBPtr;

// ----- definition for class `SphericalDSBeamformer' -----
//

/**
   @class SphericalDSBeamformer
   @usage
   1) construct this object, mb = SphericalDSBeamformer(...)
   2) set the radious of the spherical array mb.set_array_geometry(..) or
   3) set the look direction mb.set_look_direction()
   4) process each block mb.next()
   @note this implementation is based on Boaz Rafaely's letter,
   "Phase-Mode versus Delay-and-Sum Spherical Microphone Array", IEEE Signal Processing Letters, vol. 12, Oct. 2005.
*/
class SphericalDSBeamformer : public EigenBeamformer {
public:
  SphericalDSBeamformer( unsigned sampleRate, unsigned fftLen = 512, bool halfBandShift = false, unsigned NC=1, unsigned maxOrder=4, bool normalizeWeight=false, const String& nm = "SphericalDSBeamformer");
  ~SphericalDSBeamformer();
  virtual gsl_vector *calc_wng();

#ifdef ENABLE_LEGACY_BTK_API
  virtual gsl_vector *calcWNG(){ return calc_wng(); }
#endif

protected:
  virtual void calc_weights_( unsigned fbinX, gsl_vector_complex *weights );
  virtual bool calc_spherical_harmonics_at_each_position_( gsl_vector *theta_s, gsl_vector *phi_s );
};

typedef Inherit<SphericalDSBeamformer, EigenBeamformerPtr> SphericalDSBeamformerPtr;

// ----- definition for class `DualSphericalDSBeamformer' -----
//

/**
   @class DualSphericalDSBeamformer
   @usage
   1) construct this object, mb = SphericalDSBeamformer(...)
   2) set the radious of the spherical array mb.set_array_geometry(..) or
   3) set the look direction mb.set_look_direction()
   4) process each block mb.next()
   @note In addition to  SphericalDSBeamformer, this class has an object of the *normal* D&S beamformer
*/
class DualSphericalDSBeamformer : public SphericalDSBeamformer {
public:
  DualSphericalDSBeamformer( unsigned sampleRate, unsigned fftLen = 512, bool halfBandShift = false, unsigned NC=1, unsigned maxOrder=4, bool normalizeWeight=false, const String& nm = "SphericalDSBeamformer");
  ~DualSphericalDSBeamformer();

  virtual SnapShotArrayPtr snapshot_array() const { return snapshot_array_; }
  virtual BeamformerWeights* beamformer_weight_object(unsigned srcX=0) const {
    return bfweight_vec2_[srcX];
  }

#ifdef ENABLE_LEGACY_BTK_API
  virtual SnapShotArrayPtr getSnapShotArray(){ return snapshot_array(); }
#endif

protected:
  virtual void calc_weights_( unsigned fbinX, gsl_vector_complex *weights );
  virtual bool alloc_steering_unit_( int unitN=1 );

  vector<BeamformerWeights *>                   bfweight_vec2_; // weights of a normal D&S beamformer.
};

typedef Inherit<DualSphericalDSBeamformer, SphericalDSBeamformerPtr> DualSphericalDSBeamformerPtr;

// ----- definition for class DOAEstimatorSRPPSphDSB' -----
//

class DOAEstimatorSRPSphDSB : public DOAEstimatorSRPBase, public SphericalDSBeamformer {
public:
  DOAEstimatorSRPSphDSB( unsigned nBest, unsigned sampleRate, unsigned fftLen = 512, bool halfBandShift = false, unsigned NC=1, unsigned maxOrder=4, bool normalizeWeight=false, const String& nm = "DOAEstimatorSRPPSphDSB" );
  ~DOAEstimatorSRPSphDSB();
  const gsl_vector_complex* next(int frame_no = -5);
  void reset();

protected:
  virtual void  calc_steering_unit_table_();
  virtual float calc_response_power_( unsigned uttX );
};

typedef Inherit<DOAEstimatorSRPSphDSB, SphericalDSBeamformerPtr> DOAEstimatorSRPSphDSBPtr;

// ----- definition for class `SphericalHWNCBeamformer' -----
//

/**
   @class SphericalHWNCBeamformer
   @usage
   1) construct this object, mb = SphericalDSBeamformer(...)
   2) set the radious of the spherical array mb.set_array_geometry(..) or
   3) set the look direction mb.set_look_direction()
   4) process each block mb.next()
*/
class SphericalHWNCBeamformer : public EigenBeamformer {
public:
  SphericalHWNCBeamformer( unsigned sampleRate, unsigned fftLen = 512, bool halfBandShift = false, unsigned NC=1, unsigned maxOrder=4, bool normalizeWeight=false, float ratio=1.0, const String& nm = "SphericalHWNCBeamformer");
  ~SphericalHWNCBeamformer();

  virtual gsl_vector *calc_wng();
  virtual void set_wng( double ratio){ ratio_=ratio; calc_wng();}

#ifdef ENABLE_LEGACY_BTK_API
  gsl_vector *calcWNG(){ return calc_wng(); }
  void setWNG( double ratio){ set_wng(ratio);}
#endif

protected:
  virtual void calc_weights_( unsigned fbinX, gsl_vector_complex *weights );

protected:
  float ratio_;
};

typedef Inherit<SphericalHWNCBeamformer, EigenBeamformerPtr> SphericalHWNCBeamformerPtr;

// ----- definition for class `SphericalGSCBeamformer' -----
//

/**
   @class SphericalGSCBeamformer
   @usage
   1) construct this object, mb = SphericalDSBeamformer(...)
   2) set the radious of the spherical array mb.set_array_geometry(..) or
   3) set the look direction mb.set_look_direction()
   4) process each block mb.next()
   @note this implementation is based on Boaz Rafaely's letter,
   "Phase-Mode versus Delay-and-Sum Spherical Microphone Array", IEEE Signal Processing Letters, vol. 12, Oct. 2005.
*/
class SphericalGSCBeamformer : public SphericalDSBeamformer {
public:
  SphericalGSCBeamformer( unsigned sampleRate, unsigned fftLen = 512, bool halfBandShift = false, unsigned NC=1, unsigned maxOrder=4, bool normalizeWeight=false, const String& nm = "SphericalGSCBeamformer");
  ~SphericalGSCBeamformer();

  virtual const gsl_vector_complex* next(int frame_no = -5);
  virtual void reset();

  void set_look_direction(double theta, double phi);
  void set_active_weights_f(unsigned fbinX, const gsl_vector* packedWeight);

#ifdef ENABLE_LEGACY_BTK_API
  void setLookDirection(double theta, double phi){ set_look_direction(theta, phi); }
  void setActiveWeights_f( unsigned fbinX, const gsl_vector* packedWeight ){ set_active_weights_f(fbinX, packedWeight); }
#endif
};

typedef Inherit<SphericalGSCBeamformer, SphericalDSBeamformerPtr> SphericalGSCBeamformerPtr;

// ----- definition for class `SphericalHWNCGSCBeamformer' -----
//

/**
   @class SphericalHWNCGSCBeamformer
   @usage
   1) construct this object, mb = SphericalHWNCGSCBeamformer(...)
   2) set the radious of the spherical array mb.set_array_geometry(..) or
   3) set the look direction mb.set_look_direction()
   4) process each block mb.next()
*/
class SphericalHWNCGSCBeamformer : public SphericalHWNCBeamformer {
public:
  SphericalHWNCGSCBeamformer( unsigned sampleRate, unsigned fftLen = 512, bool halfBandShift = false, unsigned NC=1, unsigned maxOrder=4, bool normalizeWeight=false, float ratio=1.0, const String& nm = "SphericalHWNCGSCBeamformer");
  ~SphericalHWNCGSCBeamformer();

  virtual const gsl_vector_complex* next(int frame_no = -5);
  virtual void reset();

  void set_look_direction(double theta, double phi);
  void set_active_weights_f(unsigned fbinX, const gsl_vector* packedWeight);

#ifdef ENABLE_LEGACY_BTK_API
  void setLookDirection(double theta, double phi){ set_look_direction(theta, phi); }
  void setActiveWeights_f( unsigned fbinX, const gsl_vector* packedWeight ){ set_active_weights_f(fbinX, packedWeight); }
#endif
};

typedef Inherit<SphericalHWNCGSCBeamformer, SphericalHWNCBeamformerPtr> SphericalHWNCGSCBeamformerPtr;

// ----- definition for class `DualSphericalGSCBeamformer' -----
//

/**
   @class DualSphericalGSCBeamformer
   @usage
   1) construct this object, mb = DualSphericalGSCBeamformer(...)
   2) set the radious of the spherical array mb.set_array_geometry(..) or
   3) set the look direction mb.set_look_direction()
   4) process each block mb.next()
   @note In addition to DualSphericalGSCBeamformer, this class has an object of the *normal* D&S beamformer
*/
class DualSphericalGSCBeamformer : public SphericalGSCBeamformer {
public:
  DualSphericalGSCBeamformer( unsigned sampleRate, unsigned fftLen = 512, bool halfBandShift = false, unsigned NC=1, unsigned maxOrder=4, bool normalizeWeight=false, const String& nm = "DualSphericalGSCBeamformer");
  ~DualSphericalGSCBeamformer();

  virtual SnapShotArrayPtr snapshot_array() const {return(snapshot_array_);}
  virtual BeamformerWeights* beamformer_weight_object(unsigned srcX=0) const {
    return bfweight_vec2_[srcX];
  }

#ifdef ENABLE_LEGACY_BTK_API
  virtual SnapShotArrayPtr getSnapShotArray(){return(snapshot_array_);}
#endif

protected:
  virtual void calc_weights_( unsigned fbinX, gsl_vector_complex *weights );
  virtual bool alloc_steering_unit_( int unitN=1 );

  vector<BeamformerWeights *>                   bfweight_vec2_; // weights of a normal D&S beamformer.
};

typedef Inherit<DualSphericalGSCBeamformer, SphericalGSCBeamformerPtr> DualSphericalGSCBeamformerPtr;

// ----- definition for class `SphericalMOENBeamformer' -----
//

/**
   @class SphericalGSCBeamformer
   @usage
   1) construct this object, mb = SphericalDSBeamformer(...)
   2) set the radious of the spherical array mb.set_array_geometry(..) or
   3) set the look direction mb.set_look_direction()
   4) process each block mb.next()
   @note this implementation is based on Z. Li and R. Duraiswami's letter,
   "Flexible and Optimal Design of Spherical Microphone Arrays for Beamforming", IEEE Trans. SAP.
*/
class SphericalMOENBeamformer : public SphericalDSBeamformer {
public:
  SphericalMOENBeamformer( unsigned sampleRate, unsigned fftLen = 512, bool halfBandShift = false, unsigned NC=1, unsigned maxOrder=4, bool normalizeWeight=false, const String& nm = "SphericalMOENBeamformer");
  ~SphericalMOENBeamformer();

  virtual const gsl_vector_complex* next(int frame_no = -5);
  virtual void reset();
  void fix_terms(bool flag){ is_term_fixed_ = flag; }
  void set_diagonal_looading(unsigned fbinX, float diagonalWeight);
  virtual SnapShotArrayPtr snapshot_array() const { return snapshot_array_; }
  virtual gsl_matrix *beampattern(unsigned fbinX, double theta = 0, double phi = 0,
                                  double minTheta=-M_PI, double maxTheta=M_PI,
                                  double minPhi=-M_PI, double maxPhi=M_PI,
                                  double widthTheta=0.1, double widthPhi=0.1 );

#ifdef ENABLE_LEGACY_BTK_API
  void fixTerms( bool flag ){ fix_terms(flag); }
  void setLevelOfDiagonalLoading( unsigned fbinX, float diagonalWeight){ set_diagonal_looading(fbinX, diagonalWeight); }
  virtual gsl_matrix *getBeamPattern( unsigned fbinX, double theta = 0, double phi = 0,
                                      double minTheta=-M_PI, double maxTheta=M_PI,
                                      double minPhi=-M_PI, double maxPhi=M_PI,
                                      double widthTheta=0.1, double widthPhi=0.1 ){
    return beampattern(fbinX, theta, phi, minTheta, maxTheta, minPhi, maxPhi, widthTheta, widthPhi);
  }
#endif

protected:
  virtual bool alloc_steering_unit_( int unitN=1 );
  virtual void calc_weights_( unsigned fbinX, gsl_vector_complex *weights );
  bool calc_moen_weights_( unsigned fbinX, gsl_vector_complex *weights, double dThreshold = 1.0E-8, bool calcInverseMatrix = true, unsigned unitX=0 );

private:
  // maxOrder_ == Neff in the Li's paper.
  unsigned             bf_order_; // N in the Li's paper.
  float                CN_;
  gsl_matrix_complex** A_; /* A_[fftLen2+1][dim_][nChan]; Coeffcients of the spherical harmonics expansion; See Eq. (31) & (32) */
  gsl_matrix_complex** fixedW_; /* _fixedW[fftLen2+1][nChan][dim_]; [ A^H A + l^2 I ]^{-1} A^H */
  gsl_vector_complex** BN_;     // _BN[fftLen2+1][dim_]
  float*               diagonal_weights_;
  bool                 is_term_fixed_;
  float                dthreshold_;
};

typedef Inherit<SphericalMOENBeamformer, SphericalDSBeamformerPtr> SphericalMOENBeamformerPtr;

// ----- definition for class `SphericalSpatialDSBeamformer' -----
//

/**
   @class SphericalSpatialDSBeamformer
   @usage
   1) construct this object, mb = SphericalSpatialDSBeamformer(...)
   2) set the radious of the spherical array mb.set_array_geometry(..) or
   3) set the look direction mb.set_look_direction()
   4) process each block mb.next()
   @note this implementation is based on Boaz Rafaely's letter,
   "Phase-Mode versus Delay-and-Sum Spherical Microphone Array", IEEE Signal Processing Letters, vol. 12, Oct. 2005.
*/
class SphericalSpatialDSBeamformer : public SphericalDSBeamformer {
public:
  SphericalSpatialDSBeamformer( unsigned sampleRate, unsigned fftLen = 512, bool halfBandShift = false, unsigned NC=1, unsigned maxOrder=4, bool normalizeWeight=false, const String& nm = "SphericalSpatialDSBeamformer");
  ~SphericalSpatialDSBeamformer();
  virtual const gsl_vector_complex* next(int frame_no = -5);

protected:
  virtual void calc_weights_( unsigned fbinX, gsl_vector_complex *weights );
  virtual bool alloc_steering_unit_( int unitN = 1 );
  virtual bool calc_steering_unit_( int unitX = 0, bool isGSC = false );
};

typedef Inherit<SphericalSpatialDSBeamformer, SphericalDSBeamformerPtr> SphericalSpatialDSBeamformerPtr;

// ----- definition for class `SphericalSpatialHWNCBeamformer' -----
//

/**
   @class SphericalSpatialHWNCBeamformer
   @usage
   1) construct this object, mb = SphericalDSBeamformer(...)
   2) set the radious of the spherical array mb.set_array_geometry(..) or
   3) set the look direction mb.set_look_direction()
   4) process each block mb.next()
*/
class SphericalSpatialHWNCBeamformer : public SphericalHWNCBeamformer {
public:
  SphericalSpatialHWNCBeamformer( unsigned sampleRate, unsigned fftLen = 512, bool halfBandShift = false, unsigned NC=1, unsigned maxOrder=4, bool normalizeWeight=false, float ratio=1.0, const String& nm = "SphericalHWNCBeamformer");
  ~SphericalSpatialHWNCBeamformer();
  virtual const gsl_vector_complex* next(int frame_no = -5);

protected:
  virtual void calc_weights_( unsigned fbinX, gsl_vector_complex *weights );
  virtual bool alloc_steering_unit_( int unitN = 1 );
  virtual bool calc_steering_unit_( int unitX = 0, bool isGSC = false );

private:
  gsl_matrix_complex *calc_diffuse_noise_model_( unsigned fbinX );

  gsl_matrix_complex **SigmaSI_; // SigmaSI_[fftLen/2+1][chanN]
  double  dthreshold_;
};

typedef Inherit<SphericalSpatialHWNCBeamformer, SphericalHWNCBeamformerPtr> SphericalSpatialHWNCBeamformerPtr;

#endif


text:
/**
 * Created: Wed Aug 31 20:27:45 CST 2016
 * @author: Lie Yan
 * @email: robin.lie.yan@outlook.com
 *
 */

#pragma once

#include <gsl.h>

namespace formulae {
namespace interm {

enum class ItemTag : uint8_t { INoad, IPen, ISpace, IStyle };

struct Item {
  virtual ItemTag tag() const = 0;
  virtual gsl::owner<Item *> clone() const = 0;
};

class ItemRef {
public:
  ItemRef(gsl::owner<Item *> itemPtr) : _itemPtr(itemPtr) {}
  ItemRef(const ItemRef &rhs) { _itemPtr.reset(rhs._itemPtr->clone()); }
  ItemRef &operator=(const ItemRef &rhs) {
    _itemPtr.reset(rhs._itemPtr->clone());
    return *this;
  }
  ItemRef(ItemRef &&rhs) = default;
  ItemRef &operator=(ItemRef &&rhs) = default;
  virtual ~ItemRef() = default;

  /* -------------------------------------------------------- */
  ItemTag tag() const { return _itemPtr->tag(); }

  template<typename T>
  T& as() const { return *dynamic_cast<T*>(_itemPtr.get()); }

private:
  std::unique_ptr<Item> _itemPtr;
};

using ilist_t = list_t<ItemRef>;

}
}


text:
// Copyright @ 2012 Zhi Qiu and Chun Shen
// Ver 2.1
// Note the version used in iSS is slightly slower than the one
// used in iS, do not mix them unless necessary.

#ifndef SRC_EMISSIONFUNCTION_H_
#define SRC_EMISSIONFUNCTION_H_

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

#include <string>
#include <vector>
#include <array>
#include <memory>

#include "Table.h"
#include "TableFunction.h"
#include "ParameterReader.h"
#include "particle_decay.h"
#include "data_struct.h"
#include "data_struct.h"
#include "Random.h"
#include "pretty_ostream.h"

class EmissionFunctionArray {
 private:
    int hydro_mode;   // switch for (2+1)-d or (3+1)-d hypersurface
    int flag_PCE_;
    int flag_restrict_deltaf;
    double deltaf_max_ratio;

    pretty_ostream messager;

    const std::string path_;
    const std::string table_path_;
    const AfterburnerType afterburner_type_;

    std::shared_ptr<RandomUtil::Random> ran_gen_ptr;

    int MC_sampling;
    int number_of_repeated_sampling;
    int local_charge_conservation;

    Table *pT_tab, *phi_tab, *y_minus_eta_tab;
    int pT_tab_length, phi_tab_length;

    // the y_minus_eta_min_index holds the index
    // to the smallest positive y-eta_s value
    int y_minus_eta_tab_length, y_minus_eta_min_index;
    // true if y_minus_eta_tab has only positive part
    bool positive_y_minus_eta_table_only;

    ParameterReader *paraRdr; // used to pass-in parameters
    int USE_OSCAR_FORMAT;
    int USE_GZIP_FORMAT;
    int INCLUDE_DELTAF, INCLUDE_BULK_DELTAF, INCLUDE_DIFFUSION_DELTAF;
    int bulk_deltaf_kind;

    int turn_on_rhob;

    Table *dN_pTdpTdphidy;  // dN / (pt dpt dphi dy)
    // store the largest element when summing over xt and eta
    // to get dN / (pt dpt dphi dy); used during sampling.
    Table *dN_pTdpTdphidy_max;

    double **dN_dxtdetady;  // dN / (d^2x_t deta dy) (correction: and dtau)
    // dN / (d^2x_t deta dy) is the (weighted) sum of all elements of
    // the [pT_tab_length][phi_tab_length]-sized dN/dall matrix;
    // this array records the maximum value of this dN/dall matrix,
    // for each eta and FO-cell.
    double **dN_dxtdetady_pT_max;
    double **dN_dxtdy_4all;  // dN / (dxt dy) for all particles

    // dN/(dxt dy) for one particle species
    std::vector<double> dN_dxtdy_for_one_particle_species;

    int number_of_chosen_particles;
    // used in spectra and flow calculations;
    // it has length Nparticle, 0 means miss, 1 means include
    std::vector<int> chosen_particles_01_table;

    // 0/1: Particles with similar mass and chemical potentials
    // will be sampled using the same dN/(dxt deta dy) matrix
    int grouping_particles;

    // Usable only when grouping_particles is 1.
    // If two particles have mass and chemical potentials
    // close within this relative tolerance, they are considered to be
    // identical and will be sampled successively without regenerating
    // the dN / (dxt deta dy) matrix for efficiency.
    double mc_grouping_tolerance;

    // store particle index;
    // the sampling process follows the order specified by this table
    std::vector<int> chosen_particles_sampling_table;

    // store chosen particle monte carlo number that are not found in pdg.dat
    std::vector<int> unidentifiedPid_table;

    // list for information for all particles
    int Nparticles;
    std::vector<particle_info> particles;

    // list for information for all fluid cells
    long FO_length;
    const std::vector<FO_surf> &FOsurf_ptr;

    // store the last particle index being used by calculate_dNArrays function
    int last_particle_idx;

    double **trig_phi_table, **hypertrig_y_minus_eta_table;
    inline long determine_number_to_sample(
                    double dN, int model=1, double para1=0);

    const gsl_rng_type *gsl_type_random_number;
    gsl_rng *gsl_random_r;

    bool particles_are_the_same(int, int);
    //long *sorted_FZ;

    //array for bulk delta f coefficients
    Table *bulkdf_coeff;

    // table parameter for diffusion deltaf coefficient
    int deltaf_qmu_coeff_table_length_T;
    int deltaf_qmu_coeff_table_length_mu;
    double delta_qmu_coeff_table_T0, delta_qmu_coeff_table_mu0;
    double delta_qmu_coeff_table_dT, delta_qmu_coeff_table_dmu;
    double **deltaf_qmu_coeff_tb;       // kappa_B [1/fm^3]

    // arrays to speed up computing particle yield
    double sf_dx, sf_x_min, sf_x_max;
    int sf_tb_length;
    int sf_expint_truncate_order;
    double** sf_bessel_Kn;
    double** sf_expint_En;

    double lambert_x_min, lambert_x_max, lambert_dx;
    std::vector<double> lambert_W;

    int flag_output_samples_into_files;
    int flag_store_samples_in_memory;

    //! particle decay
    int flag_perform_decays;
    particle_decay *decayer_ptr;

 public:
    EmissionFunctionArray(std::shared_ptr<RandomUtil::Random> ran_gen,
                          Table* chosen_particle, Table* pt_tab_in,
                          Table* phi_tab_in, Table* eta_tab_in,
                          std::vector<particle_info> particles_in,
                          const std::vector<FO_surf> &FOsurf_ptr_in,
                          int flag_PCE_in, ParameterReader* paraRdr_in,
                          std::string path_in, std::string table_path,
                          AfterburnerType afterburner_type);
    ~EmissionFunctionArray();

    std::vector< std::vector<iSS_Hadron>* >* Hadron_list;

    void initialize_special_function_arrays();
    double get_special_function_K1(double arg);
    double get_special_function_K2(double arg);
    void get_special_function_En(double arg, std::vector<double> &results);
    double get_special_function_lambertW(double arg);

    void calculate_dNArrays(int);
    void calculate_dN_dxtdetady(int);
    void calculate_dN_pTdpTdphidy(int);
    void write_dN_pTdpTdphidy_toFile();
    std::string dN_pTdpTdphidy_filename;  // where to save
    void write_dN_dxtdetady_toFile();
    std::string dN_dxtdetady_filename;

    void calculate_flows(int to_order, std::string flow_diff_filename,
                         std::string flow_inte_filename);
    std::string flow_differential_filename_old, flow_integrated_filename_old;
    std::string flow_differential_filename, flow_integrated_filename;

    void calculate_dN_pTdpTdphidy_and_flows_4all_old_output(
                                                    int perform_sampling = 0);
    void calculate_dN_pTdpTdphidy_and_flows_4all(int perform_sampling = 0);
    void calculate_dN_dxtdetady_and_sample_4all();

    void shell();  // it all starts here...

    void combine_samples_to_OSCAR();
    std::string OSCAR_header_filename, OSCAR_output_filename;
    void combine_samples_to_gzip_file();

    // Sample files
    // where samples, its control informations, and its "format file"
    // are stored;
    // the first two can contain a "%d" string to generate multiple files
    std::string samples_format_filename;

    // First sampling method
    void sample_using_dN_dxtdetady_smooth_pT_phi();
    void calculate_dN_dtau_using_dN_dxtdetady(
                    double tau0 = 0, double dtau = 0.5, double tau_max = 17);
    void calculate_dN_dphi_using_dN_pTdpTdphidy();
    void calculate_dN_deta_using_dN_dxtdetady();
    void calculate_dN_dxt_using_dN_dxtdetady();
    void calculate_dN_dx_using_dN_dxtdetady(
                    double x_min, double x_max, double dx);

    // Second sampling method
    void calculate_dN_analytic(const particle_info* particle, double mu,
                               double Temperature,
                               std::array<double, 5> &results);

    // the following variables need to be set first in order to
    // call this function
    void calculate_dN_dxtdy_4all_particles();

    void calculate_dN_dxtdy_for_one_particle_species(const int particle_idx);

    double calculate_total_FZ_energy_flux();
    // to be used after calculate_dN_dxtdy_4all_particles
    void sample_using_dN_dxtdy_4all_particles_conventional();

    // Third sampling method
    void sample_using_dN_pTdpTdphidy();
    Table pT_tab4Sampling, phi_tab4Sampling;
    int pT_tab4Sampling_length, phi_tab4Sampling_length;
    double** trig_phi_tab4Sampling;

    void getbulkvisCoefficients(double Tdec,
                                std::array<double, 3> &bulkvisCoefficients);
    void load_deltaf_qmu_coeff_table(std::string filename);
    double get_deltaf_qmu_coeff(double T, double muB);

    void check_samples_in_memory();
    int get_number_of_sampled_events() {return(Hadron_list->size());};
    int get_number_of_particles(int iev) {
        return((*Hadron_list)[iev]->size());
    };

    iSS_Hadron get_hadron(int iev, int ipart) {
        return((*(*Hadron_list)[iev])[ipart]);
    };

    std::vector<iSS_Hadron>* get_hadron_list_iev(const int iev) {
        return(Hadron_list->at(iev));
    }

    void perform_resonance_feed_down(
                std::vector< std::vector<iSS_Hadron>* >* input_particle_list);
    int compute_number_of_sampling_needed(int number_of_particles_needed);
    double estimate_ideal_maximum(
        int sign, double mass, double Tdec, double mu, double f0_mass,
        TableFunction &z_exp_m_z);
    double estimate_shear_viscous_maximum(
        int sign, double mass, double Tdec, double mu, double f0_mass,
        TableFunction &z_exp_m_z, double pi_size);
    double estimate_diffusion_maximum(
        int sign, int baryon, double mass, double Tdec, double mu,
        double f0_mass, TableFunction &z_exp_m_z,
        double prefactor_qmu, double guess_ideal, double q_size);
    double get_deltaf_bulk(
        double mass, double pdotu, double bulkPi, double Tdec, int sign,
        double f0, const std::array<double, 3> bulkvisCoefficients);
    int sample_momemtum_from_a_fluid_cell(
        const double mass, const double degen, const int sign,
        const int baryon, const int strange, const int charge,
        const double pT_to, const double y_minus_eta_s_range,
        const double maximum_guess, const FO_surf *surf,
        const std::array<double, 3> bulkvisCoefficients,
        const double deltaf_qmu_coeff,
        double &pT, double &phi, double &y_minus_eta_s);
    double estimate_maximum(
        const FO_surf *surf, const int real_particle_idx, const double mass,
        const double sign, const double degen,
        const int baryon, const int strange, const int charge,
        TableFunction &z_exp_m_z,
        const std::array<double, 3> bulkvisCoefficients,
        const double deltaf_qmu_coeff);
    std::string add_one_sampled_particle(
        const int repeated_sampling_idx,
        const unsigned long FO_idx, const FO_surf *surf,
        const int particle_monval, const double mass,
        const double pT, const double phi,
        const double y_minus_eta_s, const double eta_s);
};

#endif  // SRC_EMISSIONFUNCTION_H_

/*----------------------------------------------------------------------
Change log: See changelog.txt.
----------------------------------------------------------------------*/

max_stars_repo_name: doliinychenko/iSS

max_issues_repo_name: doliinychenko/iSS

max_forks_repo_name: doliinychenko/iSS


text:
// Copyright (c) Microsoft Corporation.
// Licensed under the MIT License.

//------------------------------------------------------------------------------
// Intrinsics.h
//------------------------------------------------------------------------------

#pragma once

#include <arcana/analysis/introspector.h>

#include <gsl/gsl>

#include <array>

namespace mage
{
    class Intrinsics
    {
    public:
        Intrinsics() = default;
        Intrinsics(gsl::span<const float, 4> intrinsicsCxCyFxFy, uint32_t imageWidth, uint32_t imageHeight);
        Intrinsics(float cx, float cy, float fx, float fy, uint32_t imageWidth, uint32_t imageHeight);

        // cx, cy, fx, fy
        gsl::span<const float, 4> GetCoefficients() const { return m_coefficients; }
        std::array<float, 4> GetNormalizedCoefficients() const;

        float GetCx() const { return m_coefficients[0]; } //principle point pixels
        float GetCy() const { return m_coefficients[1]; } //principle point pixels
        float GetFx() const { return m_coefficients[2]; } //focal length pixels
        float GetFy() const { return m_coefficients[3]; } //focal length pixels

        void SetCx(float cx) { m_coefficients[0] = cx; }  //principle point pixels
        void SetCy(float cy) { m_coefficients[1] = cy; }  //principle point pixels
        void SetFx(float fx) { m_coefficients[2] = fx; }  //focal length pixels
        void SetFy(float fy) { m_coefficients[3] = fy; }  //focal length pixels

        uint32_t GetCalibrationWidth() const { return m_widthPixels; }
        uint32_t GetCalibrationHeight() const { return m_heightPixels; }

    private:
        std::array<float, 4> m_coefficients{}; // cx, cy, fx, fy in pixels
        uint32_t m_widthPixels{};  //widht/height calibration was performed for
        uint32_t m_heightPixels{};
    };

    template<typename ArchiveT>
    void introspect_object(mira::introspector<ArchiveT>& intro, const Intrinsics& intrin)
    {
        intro(
            cereal::make_nvp("Coefficients", intrin.GetCoefficients())
        );
    }
}

max_stars_repo_name: syntheticmagus/mageslam

max_issues_repo_name: syntheticmagus/mageslam

max_forks_repo_name: syntheticmagus/mageslam


text:
#include <stdio.h>
#include <stdarg.h>
#include <string.h>
#include <math.h>
#include <gbpLib.h>
#include <gbpMisc.h>
#include <gbpRNG.h>
#include <gsl/gsl_linalg.h>

void add_gaussian_noise(double *data, int n_data, int *seed, double sigma, double *covariance) {
    int       i_data;
    int       j_data;
    double *  data_temp;
    RNG_info *RNG;

    // Initialize RNG
    RNG = (RNG_info *)SID_malloc(sizeof(RNG_info));
    init_RNG(seed, RNG, RNG_DEFAULT);

    // Generate a random displacement vector
    gsl_vector *b;
    b = gsl_vector_calloc(n_data);
    for(i_data = 0; i_data < n_data; i_data++)
        gsl_vector_set(b, i_data, sigma * random_gaussian(RNG));

    // Use the rotated covariance matrix (if available)
    if(covariance != NULL) {
        // Perform Cholesky Decomposition
        gsl_matrix *m;
        m = gsl_matrix_calloc(n_data, n_data);
        for(i_data = 0; i_data < n_data; i_data++)
            for(j_data = 0; j_data < n_data; j_data++)
                gsl_matrix_set(m, i_data, j_data, covariance[i_data * n_data + j_data]);
        gsl_linalg_cholesky_decomp(m);
        for(i_data = 0; i_data < n_data; i_data++)
            for(j_data = i_data + 1; j_data < n_data; j_data++)
                gsl_matrix_set(m, i_data, j_data, 0.);
        for(i_data = 0; i_data < n_data; i_data++) {
            for(j_data = 0; j_data < n_data; j_data++)
                data[i_data] += gsl_matrix_get(m, i_data, j_data) * gsl_vector_get(b, j_data);
        }
        gsl_matrix_free(m);
    } else {
        for(j_data = 0; j_data < n_data; j_data++)
            data[j_data] += gsl_vector_get(b, j_data);
    }

    // Clean-up
    free_RNG(RNG);
    SID_free(SID_FARG RNG);
    gsl_vector_free(b);
}

max_stars_repo_name: gbpoole/gbpCode

max_issues_repo_name: gbpoole/gbpCode

max_forks_repo_name: gbpoole/gbpCode


text:
//
//  Engine.hpp
//  PretendToWork
//
//  Created by tqtifnypmb on 14/12/2017.
//  Copyright © 2017 tqtifnypmb. All rights reserved.
//

#pragma once

#include "../rope/Rope.h"
#include "../types.h"
#include "Revision.h"

#include <gsl/gsl>
#include <vector>
#include <map>

namespace brick
{

class Engine {
public:
    using DeltaList = std::vector<Revision>;

    Engine(size_t authorId, gsl::not_null<Rope*> rope);

    template <class Converter>
    void insert(gsl::span<const char> bytes, size_t pos);
    void insert(const detail::CodePointList& cplist, size_t pos);

    void erase(const Range& range);

    std::pair<DeltaList, DeltaList> sync(Engine& other);

    void fastForward(const std::vector<Revision>& revs);

    const std::vector<Revision>& revisions() const {
        return revisions_;
    }

    std::vector<Revision>& revisions() {
        return revisions_;
    }

    size_t authorId() const {
        return authorId_;
    }

    size_t nextRevId() {
        return revId_++;
    }

private:
    Revision delta(const Revision& history, Revision& rev);
    std::vector<Revision> delta(Revision& rev);

    void appendRevision(Revision rev);
    bool appendRevision(Revision rev, bool pendingRev, std::vector<Revision>* deltas);
    std::pair<DeltaList, DeltaList> sync(Engine& other, bool symetric);

    gsl::not_null<Rope*> rope_;
    std::vector<Revision> revisions_;
    std::vector<Revision> pendingRevs_;
    std::map<size_t, std::pair<size_t, size_t>> sync_state_;

    size_t authorId_;
    size_t revId_;
};

template <class Converter>
void Engine::insert(gsl::span<const char> bytes, size_t pos) {
    return insert(Converter::encode(bytes), pos);
}

}   // namespace brick

max_stars_repo_name: tqtifnypmb/brick

max_issues_repo_name: tqtifnypmb/brick

max_forks_repo_name: tqtifnypmb/brick


text:
/* linalg/qr_ud.c
 *
 * Copyright (C) 2020 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 <stdlib.h>
#include <string.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>

/*
 * this module contains routines for the QR factorization of a matrix
 * using the recursive Level 3 BLAS algorithm of Elmroth and Gustavson with
 * additional modifications courtesy of Julien Langou.
 */

static int aux_QR_TRD_decomp (gsl_matrix * U, gsl_matrix * A, const gsl_vector * D, gsl_matrix * T);
static double qrtd_householder_transform (double *v0, double *v1);
static double qrtrd_householder_transform (double *v0, gsl_vector * v, double *d);

/*
gsl_linalg_QR_UD_decomp()
  Compute the QR decomposition of the "triangle on top of diagonal" matrix

  [ U ] = Q [ R ]
  [ D ]     [ 0 ]

where U is N-by-N upper triangular, D is N-by-N diagonal

Inputs: U - on input, upper triangular N-by-N matrix
            on output, R factor in upper triangle
        D - diagonal N-by-N matrix
        Y - (output) upper triangular N-by-N matrix (lower part of V)
        T - (output) block reflector matrix, N-by-N

Notes:
1) Based on the Elmroth/Gustavson algorithm, taking into account the
sparse structure of the U,S matrices

2) The Householder matrix V has the special form:

      N
V = [ I ] N
    [ Y ] N

with Y upper triangular

3) The orthogonal matrix is

Q = I - V T V^T
*/

int
gsl_linalg_QR_UD_decomp (gsl_matrix * U, const gsl_vector * D, gsl_matrix * Y, gsl_matrix * T)
{
  const size_t N = U->size1;

  if (N != U->size2)
    {
      GSL_ERROR ("U matrix must be square", GSL_ENOTSQR);
    }
  else if (D->size != N)
    {
      GSL_ERROR ("D matrix does not match dimensions of U", GSL_EBADLEN);
    }
  else if (Y->size1 != Y->size2)
    {
      GSL_ERROR ("Y matrix must be square", GSL_ENOTSQR);
    }
  else if (Y->size1 != N)
    {
      GSL_ERROR ("Y matrix does not match dimensions of U", GSL_EBADLEN);
    }
  else if (T->size1 != N || T->size2 != N)
    {
      GSL_ERROR ("T matrix has wrong dimensions", GSL_EBADLEN);
    }
  else if (N == 1)
    {
      /* base case, compute Householder transform for single column matrix */
      double * T00 = gsl_matrix_ptr(T, 0, 0);
      double * U00 = gsl_matrix_ptr(U, 0, 0);
      double * Y00 = gsl_matrix_ptr(Y, 0, 0);
      *Y00 = gsl_vector_get(D, 0);
      *T00 = qrtd_householder_transform(U00, Y00);
      return GSL_SUCCESS;
    }
  else
    {
      /*
       * partition matrices:
       *
       *       N1  N2              N1  N2
       * N1 [ U11 U12 ] and  N1 [ T11 T12 ]
       * N2 [  0  U22 ]      N2 [  0  T22 ]
       * N1 [ D11 D12 ]
       * N2 [  0  D22 ]
       */
      int status;
      const size_t N1 = N / 2;
      const size_t N2 = N - N1;

      gsl_matrix_view U11 = gsl_matrix_submatrix(U, 0, 0, N1, N1);
      gsl_matrix_view U12 = gsl_matrix_submatrix(U, 0, N1, N1, N2);
      gsl_matrix_view U22 = gsl_matrix_submatrix(U, N1, N1, N2, N2);

      gsl_vector_const_view D1 = gsl_vector_const_subvector(D, 0, N1);
      gsl_vector_const_view D2 = gsl_vector_const_subvector(D, N1, N2);

      gsl_matrix_view Y11 = gsl_matrix_submatrix(Y, 0, 0, N1, N1);
      gsl_matrix_view Y12 = gsl_matrix_submatrix(Y, 0, N1, N1, N2);

      gsl_matrix_view T11 = gsl_matrix_submatrix(T, 0, 0, N1, N1);
      gsl_matrix_view T12 = gsl_matrix_submatrix(T, 0, N1, N1, N2);
      gsl_matrix_view T22 = gsl_matrix_submatrix(T, N1, N1, N2, N2);

      gsl_matrix_view m;

      /*
       * Eq. 2: recursively factor
       *
       *       N1           N1
       * N1 [ U11 ] = Q1 [ R11 ] N1
       * N2 [  0  ]      [  0  ] N2
       * N1 [ D11 ]      [  0  ] N1
       * N2 [  0  ]      [  0  ] N2
       */
      status = gsl_linalg_QR_UD_decomp(&U11.matrix, &D1.vector, &Y11.matrix, &T11.matrix);
      if (status)
        return status;

      /*
       * Eq. 3:
       *
       *      N2              N2            N2
       * N1 [ R12  ] = Q1^T [ U12 ] = [   U12 - W   ] N1
       * N2 [ U22~ ]        [ U22 ]   [     U22     ] N2
       * N1 [ D12~ ]        [  0  ]   [  - V31 W    ] N1
       * N2 [ D22~ ]        [ D22 ]   [     D22     ] N2
       *
       * where W = T11^T U12, using T12 as temporary storage, and
       *
       *        N1
       * V1 = [  I  ] N1
       *      [  0  ] N2
       *      [ V31 ] N1
       *      [  0  ] N2
       */
      gsl_matrix_memcpy(&T12.matrix, &U12.matrix);                                                    /* W := U12 */
      gsl_blas_dtrmm(CblasLeft, CblasUpper, CblasTrans, CblasNonUnit, 1.0, &T11.matrix, &T12.matrix); /* W := T11^T W */
      gsl_matrix_sub(&U12.matrix, &T12.matrix);                                                       /* R12 := U12 - W */
      gsl_blas_dtrmm(CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, -1.0, &Y11.matrix, &T12.matrix); /* W := -V31 W */
      gsl_matrix_memcpy(&Y12.matrix, &T12.matrix);                                                    /* Y12 := -V31 W */

      /*
       * Eq. 4: factor the "triangle on top of rectangle on top of diagonal" matrix
       *
       * N2 [ U22 ] = Q2~ [ R22 ] N2
       * N1 [ Y12 ]       [  0  ] N1
       * N2 [ D22 ]       [  0  ] N2
       */
      m = gsl_matrix_submatrix(Y, 0, N1, N, N2);
      status = aux_QR_TRD_decomp(&U22.matrix, &m.matrix, &D2.vector, &T22.matrix);
      if (status)
        return status;

      /*
       * Eq. 13: update T12 := -T11 * V1^T * V2 * T22
       *
       * where:
       *
       *        N1                N2
       * V1 = [  I   ] N1   V2 = [  0   ] N1
       *      [  0   ] N2        [  I   ] N2
       *      [ V31~ ] N1        [ V32~ ] N1
       *      [  0   ] N2        [ V42~ ] N2
       *
       * Note: V1^T V2 = V31~^T V32~
       */

      gsl_matrix_memcpy(&T12.matrix, &Y12.matrix);                                                    /* T12 := V32~ */
      gsl_blas_dtrmm(CblasLeft, CblasUpper, CblasTrans, CblasNonUnit, 1.0, &Y11.matrix, &T12.matrix); /* T12 := V31~^T V32~ */
      gsl_blas_dtrmm(CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, -1.0, &T11.matrix, &T12.matrix); /* T12 := -T11 * T12 */
      gsl_blas_dtrmm(CblasRight, CblasUpper, CblasNoTrans, CblasNonUnit, 1.0, &T22.matrix, &T12.matrix); /* T12 := T12 * T22 */

      return GSL_SUCCESS;
    }
}

/* Find the least squares solution to the overdetermined system
 *
 *   [ U ] x = b
 *   [ D ]
 *
 * using the QR factorization [ U; D ] = Q R.
 *
 * Inputs: R    - upper triangular R matrix, N-by-N
 *         Y    - upper triangular Y matrix, N-by-N
 *         T    - upper triangular block reflector, N-by-N
 *         b    - right hand side, size 2*N
 *         x    - (output) solution, size 2*N
 *                x(1:N) = least squares solution vector
 *                x(N+1:2*N) = vector whose norm equals ||b - Ax||
 *         work - workspace, size N
 */

int
gsl_linalg_QR_UD_lssolve (const gsl_matrix * R, const gsl_matrix * Y, const gsl_matrix * T,
                          const gsl_vector * b, gsl_vector * x, gsl_vector * work)
{
  const size_t N = R->size1;
  const size_t M = 2 * N;

  if (R->size2 != N)
    {
      GSL_ERROR ("R matrix must be square", GSL_ENOTSQR);
    }
  else if (Y->size1 != Y->size2)
    {
      GSL_ERROR ("Y matrix must be square", GSL_ENOTSQR);
    }
  else if (Y->size1 != N)
    {
      GSL_ERROR ("Y and R must have same dimensions", GSL_EBADLEN);
    }
  else if (T->size1 != N || T->size2 != N)
    {
      GSL_ERROR ("T matrix must be N-by-N", GSL_EBADLEN);
    }
  else if (M != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (M != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else if (N != work->size)
    {
      GSL_ERROR ("workspace must be length N", GSL_EBADLEN);
    }
  else
    {
      return gsl_linalg_QR_UU_lssolve(R, Y, T, b, x, work);
    }
}

/*
aux_QR_TRD_decomp()
  Compute the QR decomposition of the "triangle on top of rectangle on top of diagonal" matrix

  [ U ] = Q [ R ]
  [ A ]     [ 0 ]
  [ D ]     [ 0 ]

where U is N-by-N upper triangular, A is M-by-N dense, D is N-by-N diagonal

Inputs: U - on input, upper triangular N-by-N matrix
            on output, R factor in upper triangle
        A - on input, A is stored in A(1:M,1:N)
            on output, Y is stored in A(1:M+N,1:N)
            size (M+N)-by-N
        D - diagonal N-by-N matrix
        T - (output) block reflector matrix, N-by-N

Notes:
1) Based on the Elmroth/Gustavson algorithm, taking into account the
sparse structure of the U,S matrices

2) The Householder matrix V has the special form:

      N
V = [ I ] N
    [ Y ] M+N

with Y upper trapezoidal

3) The orthogonal matrix is

Q = I - V T V^T
*/

static int
aux_QR_TRD_decomp (gsl_matrix * U, gsl_matrix * A, const gsl_vector * D, gsl_matrix * T)
{
  const size_t N = U->size1;

  if (N != U->size2)
    {
      GSL_ERROR ("U matrix must be square", GSL_ENOTSQR);
    }
  else if (A->size1 <= N)
    {
      GSL_ERROR ("A matrix must have > N rows", GSL_EBADLEN);
    }
  else if (D->size != N)
    {
      GSL_ERROR ("D matrix does not match dimensions of U", GSL_EBADLEN);
    }
  else if (T->size1 != N || T->size2 != N)
    {
      GSL_ERROR ("T matrix has wrong dimensions", GSL_EBADLEN);
    }
  else if (N == 1)
    {
      /* base case, compute Householder transform for single column matrix */
      const size_t M = A->size1 - N;
      double * T00 = gsl_matrix_ptr(T, 0, 0);
      double * U00 = gsl_matrix_ptr(U, 0, 0);
      gsl_vector_view v = gsl_matrix_subcolumn(A, 0, 0, M);
      double * D00 = gsl_matrix_ptr(A, M, 0);
      *D00 = gsl_vector_get(D, 0);
      *T00 = qrtrd_householder_transform(U00, &v.vector, D00);
      return GSL_SUCCESS;
    }
  else
    {
      /*
       * partition matrices:
       *
       *       N1  N2              N1  N2
       * N1 [ U11 U12 ] and  N1 [ T11 T12 ]
       * N2 [  0  U22 ]      N2 [  0  T22 ]
       * N1 [ D11 D12 ]
       * N2 [  0  D22 ]
       */
      int status;
      const size_t M = A->size1 - N;
      const size_t N1 = N / 2;
      const size_t N2 = N - N1;

      gsl_matrix_view U11 = gsl_matrix_submatrix(U, 0, 0, N1, N1);
      gsl_matrix_view U12 = gsl_matrix_submatrix(U, 0, N1, N1, N2);
      gsl_matrix_view U22 = gsl_matrix_submatrix(U, N1, N1, N2, N2);

      gsl_matrix_view A1 = gsl_matrix_submatrix(A, 0, 0, M, N1);
      gsl_matrix_view A2 = gsl_matrix_submatrix(A, 0, N1, M, N2);

      gsl_vector_const_view D1 = gsl_vector_const_subvector(D, 0, N1);
      gsl_vector_const_view D2 = gsl_vector_const_subvector(D, N1, N2);

      gsl_matrix_view T11 = gsl_matrix_submatrix(T, 0, 0, N1, N1);
      gsl_matrix_view T12 = gsl_matrix_submatrix(T, 0, N1, N1, N2);
      gsl_matrix_view T22 = gsl_matrix_submatrix(T, N1, N1, N2, N2);

      gsl_matrix_view Y41 = gsl_matrix_submatrix(A, M, 0, N1, N1);
      gsl_matrix_view Y42 = gsl_matrix_submatrix(A, M, N1, N1, N2);

      gsl_matrix_view m;

      /*
       * Eq. 2: recursively factor
       *
       *       N1           N1
       * N1 [ U11 ] = Q1 [ R11 ] N1
       * N2 [  0  ]      [  0  ] N2
       * M  [ A1  ]      [  0  ] M
       * N1 [ D11 ]      [  0  ] N1
       * N2 [  0  ]      [  0  ] N2
       */
      m = gsl_matrix_submatrix(A, 0, 0, M + N1, N1);
      status = aux_QR_TRD_decomp(&U11.matrix, &m.matrix, &D1.vector, &T11.matrix);
      if (status)
        return status;

      /*
       * Eq. 3:
       *
       *      N2              N2            N2
       * N1 [ R12  ] = Q1^T [ U12 ] = [   U12 - W   ] N1
       * N2 [ U22~ ]        [ U22 ]   [     U22     ] N2
       * M  [ A2~  ]        [ A2  ]   [ A2 - V31 W  ] M
       * N1 [ D12~ ]        [  0  ]   [  - V41 W    ] N1
       * N2 [ D22~ ]        [ D22 ]   [     D22     ] N2
       *
       * where W = T11^T ( U12 + V31^T A2 ), using T12 as temporary storage, and
       *
       *        N1
       * V1 = [  I  ] N1
       *      [  0  ] N2
       *      [ V31 ] M
       *      [ V41 ] N1
       *      [  0  ] N2
       */
      gsl_matrix_memcpy(&T12.matrix, &U12.matrix);                                                    /* W := U12 */
      gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, &A1.matrix, &A2.matrix, 1.0, &T12.matrix);        /* W := U12 + V31^T A2 */
      gsl_blas_dtrmm(CblasLeft, CblasUpper, CblasTrans, CblasNonUnit, 1.0, &T11.matrix, &T12.matrix); /* W := T11^T W */
      gsl_matrix_sub(&U12.matrix, &T12.matrix);                                                       /* R12 := U12 - W */
      gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, -1.0, &A1.matrix, &T12.matrix, 1.0, &A2.matrix);     /* A2 := A2 - V31 W */

      gsl_blas_dtrmm(CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, -1.0, &Y41.matrix, &T12.matrix); /* W := -V41 W */
      gsl_matrix_memcpy(&Y42.matrix, &T12.matrix);                                                       /* V42 := -431 W */

      /*
       * Eq. 4: factor the "triangle on top of rectangle on top of diagonal" matrix
       *
       * N2 [ U22 ] = Q2~ [ R22 ] N2
       * N1 [ Y12 ]       [  0  ] N1
       * N2 [ D22 ]       [  0  ] N2
       */
      m = gsl_matrix_submatrix(A, 0, N1, M + N, N2);
      status = aux_QR_TRD_decomp(&U22.matrix, &m.matrix, &D2.vector, &T22.matrix);
      if (status)
        return status;

      /*
       * Eq. 13: update T12 := -T11 * V1^T * V2 * T22
       *
       * where:
       *
       *        N1                N2
       * V1 = [  I   ] N1   V2 = [  0   ] N1
       *      [  0   ] N2        [  I   ] N2
       *      [ V31~ ] M         [ V32~ ] M
       *      [ V41~ ] N1        [ V42~ ] N1
       *      [  0   ] N2        [ V52~ ] N2
       *
       * Note: V1^T V2 = V31~^T V32~ + V41~^T V42~
       */

      gsl_matrix_memcpy(&T12.matrix, &Y42.matrix);                                                    /* T12 := V42~ */
      gsl_blas_dtrmm(CblasLeft, CblasUpper, CblasTrans, CblasNonUnit, 1.0, &Y41.matrix, &T12.matrix); /* T12 := V41~^T V42~ */
      gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, &A1.matrix, &A2.matrix, 1.0, &T12.matrix);        /* T12 := V31~^T V32~^T + V41~^T V42~ */
      gsl_blas_dtrmm(CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, -1.0, &T11.matrix, &T12.matrix); /* T12 := -T11 * T12 */
      gsl_blas_dtrmm(CblasRight, CblasUpper, CblasNoTrans, CblasNonUnit, 1.0, &T22.matrix, &T12.matrix); /* T12 := T12 * T22 */

      return GSL_SUCCESS;
    }
}

/*
qrtd_householder_transform()
  This routine is an optimized version of
gsl_linalg_householder_transform(), designed for the QR
decomposition of M-by-N matrices of the form:

B = [ U ]
    [ S ]

where U,S are N-by-N upper triangular.
This routine computes a householder transformation (tau,v) of a
x so that P x = [ I - tau*v*v' ] x annihilates x(1:n-1). x will
be a subcolumn of the matrix B, and so its structure will be:

x = [ x0 ] <- 1 nonzero value for the diagonal element of U
    [ y0 ] <- 1 nonzero value for the diagonal element of S

Inputs: v0 - pointer to diagonal element of U
             on input, v0 = x0;
        v1 - on input, pointer to diagonal element of S
             on output, householder vector v
*/

static double
qrtd_householder_transform (double *v0, double *v1)
{
  /* replace v[0:M-1] with a householder vector (v[0:M-1]) and
     coefficient tau that annihilate v[1:M-1] */

  double alpha, beta, tau ;

  /* compute xnorm = || [ 0 ; v ] ||, ignoring zero part of vector */
  double xnorm = fabs(*v1);

  if (xnorm == 0)
    {
      return 0.0; /* tau = 0 */
    }

  alpha = *v0;
  beta = - GSL_SIGN(alpha) * hypot(alpha, xnorm) ;
  tau = (beta - alpha) / beta ;

  {
    double s = (alpha - beta);

    if (fabs(s) > GSL_DBL_MIN)
      {
        *v1 /= s;
        *v0 = beta;
      }
    else
      {
        *v1 *= GSL_DBL_EPSILON / s;
        *v1 /= GSL_DBL_EPSILON;
        *v0 = beta;
      }
  }

  return tau;
}

/*
qrtrd_householder_transform()
  This routine is an optimized version of
gsl_linalg_householder_transform(), designed for the QR
decomposition of M-by-N matrices of the form:

B = [ U ]
    [ A ]
    [ D ]

where U is N-by-N upper triangular A is M-by-N dense and D is N-by-N diagonal.
This routine computes a householder transformation (tau,v) of a
x so that P x = [ I - tau*v*v' ] x annihilates x(1:n-1). x will
be a subcolumn of the matrix B, and so its structure will be:

x = [ x0 ] <- 1 nonzero value for the diagonal element of S
    [ 0  ] <- N - j - 1 zeros, where j is column of matrix in [0,N-1]
    [ x  ] <- M nonzero values for the dense part A

Inputs: v0 - pointer to diagonal element of S
             on input, v0 = x0;
        v  - on input, x vector
             on output, householder vector v
        d  - on input, diagonal element of D
             on output, householder element
*/

static double
qrtrd_householder_transform (double *v0, gsl_vector * v, double *d)
{
  /* replace v[0:M-1] with a householder vector (v[0:M-1]) and
     coefficient tau that annihilate v[1:M-1] */

  double alpha, beta, tau ;

  /* compute xnorm = || [ 0 ; v ; 0 ; d] ||, ignoring zero part of vector */
  double xnorm;

  xnorm = gsl_blas_dnrm2(v);
  xnorm = gsl_hypot(xnorm, *d);

  if (xnorm == 0)
    {
      return 0.0; /* tau = 0 */
    }

  alpha = *v0;
  beta = - GSL_SIGN(alpha) * hypot(alpha, xnorm) ;
  tau = (beta - alpha) / beta ;

  {
    double s = (alpha - beta);

    if (fabs(s) > GSL_DBL_MIN)
      {
        gsl_blas_dscal (1.0 / s, v);
        *d /= s;
        *v0 = beta;
      }
    else
      {
        gsl_blas_dscal (GSL_DBL_EPSILON / s, v);
        gsl_blas_dscal (1.0 / GSL_DBL_EPSILON, v);
        *d *= GSL_DBL_EPSILON / s;
        *d /= GSL_DBL_EPSILON;
        *v0 = beta;
      }
  }

  return tau;
}

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


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

#pragma once

#ifndef zip_efeefa3a_157a_4c4c_83b0_99ec8e0364a1_h
#define zip_efeefa3a_157a_4c4c_83b0_99ec8e0364a1_h

#include <gslib/config.h>
#include <gslib/tree.h>
#include <gslib/vdir.h>

__gslib_begin__

typedef void* zip_handle;

struct zip_argument
{
    int             level;      /* 0-9, default -1 */
    const gchar*    name;
    const gchar*    comment;

    zip_argument() { level = -1, comment = name = 0; }
    zip_argument(const gchar* n) { level = -1, name = n, comment = 0; }
    zip_argument(int l, const gchar* n, const gchar* c) { level = l, name = n, comment = c; }
};

struct zip_node
{
    enum
    {
        zt_folder,
        zt_file,
        zt_buffer,
    };
    uint            tag;
    string          name;
};

typedef _treenode_wrapper<zip_node> zip_wrapper;
typedef tree<zip_node, zip_wrapper> zip_directory;

struct zip_folder:
    public zip_node
{
    zip_folder() { tag = zt_folder; }
};

struct zip_file:
    public zip_node
{
    string          password;
    string          comment;
    string          path;
    zip_file() { tag = zt_file; }
};

struct zip_buffer:
    public zip_node
{
    string          password;
    string          comment;
    vessel          buffer;
    zip_buffer() { tag = zt_buffer; }
};

extern void get_zip_path(string& path, const zip_wrapper* w);
extern void create_zip_directory(zip_directory& dir, const gchar* path);
extern void convert_zip_directory(zip_directory& dir, vdir& vd);
extern bool do_zip(const zip_argument& arg, const zip_directory& dir);
extern bool do_zip(const zip_argument& arg, const gchar* path);

struct zip_info_node
{
    string          name;
    int             length;
    int             size;
    float           ratio;
    uint            crc;
};

typedef list<zip_info_node> zip_info;

extern bool view_zip_info(zip_info& info, const gchar* path);
extern void do_unzip(const gchar* src, const gchar* dest);
extern void do_unzip(const gchar* src, vdir& vd);

class zip
{
public:
    zip() { _zh = 0; }
    ~zip() { close(0); }
    void create(const zip_argument& arg);
    void close(const gchar* comment);
    bool start(const gchar* name, const gchar* password, const gchar* comment);
    bool write(const void* buf, int size);
    void finish();

protected:
    zip_handle      _zh;
    zip_argument    _arg;
};

class unzip
{
public:
    unzip() { _zh = 0; }
    ~unzip() { close(); }
    bool open(const gchar* filename);
    void close();
    bool view(zip_info& info);
    bool extract(const gchar* path);
    bool extract_item(const gchar* path, const gchar* name, const gchar* password);
    bool extract(vdir& vd);
    bool extract_item(vdir& vd, const gchar* name, const gchar* password);

protected:
    zip_handle      _zh;
};

__gslib_end__

#endif


text:
/*
 Copyright (C) 2002 M. Marques, A. Castro, A. Rubio, G. Bertsch

 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, 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.

 $Id$
*/

#include <config.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>
#include <fortran_types.h>

#include <gsl/gsl_sf_legendre.h>

/* normalization constants for the spherical harmonics */
static double sph_cnsts[9] = {
	0.282094791773878, /* l = 0 */
	0.488602511902920, 0.488602511902920, 0.488602511902920, /* l = 1 */
	0.182091405098680, 0.364182810197360, 0.630783130505040, 0.364182810197360, 0.182091405098680
};

/* Computes real spherical harmonics ylm in the direction of vector r:
	 ylm = c * plm( cos(theta) ) * sin(m*phi)   for   m <  0
	 ylm = c * plm( cos(theta) ) * cos(m*phi)   for   m >= 0
	 with (theta,phi) the polar angles of r, c a positive normalization
	 constant and plm associated legendre polynomials. */

void FC_FUNC_(oct_ylm, OCT_YLM)
     (const fint * np, const double *x, const double *y, const double *z, const fint *l, const fint *m, double * restrict ylm)
{
  double r, r2, rr, rx, ry, rz, cosphi, sinphi, cosm, sinm, phase;
  int i, ip;

  if(l[0] == 0) {
    for (ip = 0; ip < np[0]; ip++) ylm[ip] = sph_cnsts[0];
    return;
  }

  for (ip = 0; ip < np[0]; ip++){

    r2 = x[ip]*x[ip] + y[ip]*y[ip] + z[ip]*z[ip];

    /* if r=0, direction is undefined => make ylm=0 except for l=0 */
    if(r2 < 1.0e-15){
      ylm[ip] = 0.0;
      continue;
    }

    switch(l[0]){
    case 1:
      rr = 1.0/sqrt(r2);
      switch(m[0]){
      case -1: ylm[ip] = -sph_cnsts[1]*rr*y[ip]; continue;
      case  0: ylm[ip] =  sph_cnsts[2]*rr*z[ip]; continue;
      case  1: ylm[ip] = -sph_cnsts[3]*rr*x[ip]; continue;
      }
    case 2:
      switch(m[0]){
      case -2: ylm[ip] =  sph_cnsts[4]*6.0*x[ip]*y[ip]/r2;
	continue;
      case -1: ylm[ip] = -sph_cnsts[5]*3.0*y[ip]*z[ip]/r2;
	continue;
      case  0: ylm[ip] =  sph_cnsts[6]*0.5*(3.0*z[ip]*z[ip]/r2 - 1.0);
	continue;
      case  1: ylm[ip] = -sph_cnsts[7]*3.0*x[ip]*z[ip]/r2;
	continue;
      case  2: ylm[ip] =  sph_cnsts[8]*3.0*(x[ip]*x[ip] - y[ip]*y[ip])/r2;
	continue;
      }
    }

    /* get phase */
    rr = 1.0/sqrt(r2);

    rx = x[ip]*rr;
    ry = y[ip]*rr;
    rz = z[ip]*rr;

    r = hypot(rx, ry);
    if(r < 1e-20) r = 1e-20; /* one never knows... */

    cosphi = rx/r;
    sinphi = ry/r;

    /* compute sin(mphi) and cos(mphi) by adding cos/sin */
    cosm = 1.; sinm=0.;
    for(i = 0; i < abs(m[0]); i++){
      double a = cosm, b = sinm;
      cosm = a*cosphi - b*sinphi;
      sinm = a*sinphi + b*cosphi;
    }
    phase = m[0] < 0 ? sinm : cosm;
    phase = m[0] == 0 ? phase : sqrt(2.0)*phase;

    /* adding small number (~= 10^-308) to avoid floating invalids */
    rz = rz + DBL_MIN;

    r = gsl_sf_legendre_sphPlm(l[0], abs(m[0]), rz);

    /* I am not sure whether we are including the Condon-Shortley factor (-1)^m */
    ylm[ip] = r*phase;
  }
}

max_stars_repo_name: neelravi/octopus

max_issues_repo_name: neelravi/octopus

max_forks_repo_name: neelravi/octopus


text:
#pragma once

#include <string>
#include <vector>
#include <set>
#include <map>

#include <gsl/gsl>

#include <gamesystems/map/sector.h>

#include "streams.h"
#include "animgoals/anim.h"

struct QuestSave {
	GameTime acceptedTime;
	uint32_t acceptedArea;
	uint32_t state;
};

struct ObjRefSave {
	ObjectId id;
	uint32_t mapId;
	locXY location;
};

struct AnimGoalSave {
	AnimGoalType type;
	ObjRefSave selfObj;
	ObjRefSave targetObj;
	ObjRefSave blockObj;
	ObjRefSave scratchObj;
	ObjRefSave parentObj;
	locXY targetLoc;
	locXY range;
	int animId;
	int animIdPrev;
	int animData;
	int animData2;
	int spellData;
	int skillData;
	int flagsData;
	int scratchVal1;
	int scratchVal2;
	int scratchVal3;
	int scratchVal4;
	int scratchVal5;
	int scratchVal6;
	int soundHandle;
};

struct AnimSlotSave {
	AnimSlotId id;
	uint32_t flags;
	int currentState;
	uint32_t field14;
	uint64_t field18;
	ObjRefSave animObj;
	std::vector<AnimGoalSave> goals;
	std::array<uint8_t, 0xF8> unk;
	uint64_t pathInfo;
	GameTime someTime;
	uint32_t field2c90;
};

struct SpellObjSave {
	ObjectId id;
	uint32_t casterPartSysName;
};

struct SpellSave {
	uint32_t spellId;
	uint32_t spellIdx;
	uint32_t active;

	uint32_t spellEnum;
	uint32_t spellEnumOrg;
	uint32_t flagsSth;
	ObjectId caster;
	uint32_t casterPartSysName; // 0 if none
	uint32_t spellClass;
	uint32_t spellLevel;
	uint32_t baseCasterLevel;
	uint32_t spellDc;
	uint32_t spellObjCount; // Guess
	objHndl unknownObj;
	std::array<SpellObjSave, 128> spellObjs; // Not really sure if this is spell objs
	size_t targetListNumItemsCopy;
	size_t targetListNumItems;
	std::array<SpellObjSave, 32> targets;
	uint32_t numProjectiles;
	std::array<ObjectId, 5> projectiles;
	LocFull aoeCenter;
	uint32_t duration;
	uint32_t durationRemaining;
	uint32_t spellRange;
	uint32_t savingThrowResult;
	uint32_t metaMagic;
	uint32_t spellId2;
};

struct TimeEventArgSave {
	int32_t intVal = 0;
	float floatVal = 0;
	ObjRefSave objectRefVal;
	locXY locVal = {0, 0};
	std::string pythonObj;
};

struct TimeEventSave {
	GameTime expiresAt;
	TimeEventType type;
	std::vector<TimeEventArgSave> args;
};

class SaveGameArchive {
public:
	static void Create(const std::string &folder,
		const std::string &filename);

	static void Unpack(const std::string &filename,
		const std::string &folder);

private:
	static void AddFolder(const std::string &folder,
		OutputStream &indexStream,
		OutputStream &dataStream);

	static void UnpackFolder(const std::string &folder,
		InputStream &indexStream,
		InputStream &dataStream);

	static void UnpackFile(const std::string &folder,
		InputStream &indexStream,
		InputStream &dataStream);
};

class SaveGame {
public:

	void Load(const std::string &path);

private:

	void LoadDataSav(InputStream &in);

	std::vector<std::string> mCustomDescriptions;
	bool mIronman = false;
	uint32_t mIronmanSaveNumber = 0;
	std::string mIronmanSaveName;
	std::vector<SectorTime> mSectorTimes;
	std::vector<uint8_t> mGlobalVars;
	std::vector<uint8_t> mGlobalFlags;
	uint32_t mStoryState = 0;
	std::vector<int32_t> mEncounterQueue;

	std::string mMapDataPath;
	std::string mMapSavePath;
	std::set<uint32_t> mVisitedMaps;

	uint32_t mNextSpellId = 0;
	std::vector<SpellSave> mSpells;

	uint32_t mLightSchemeId = 0;
	uint32_t mLightSchemeHourOfDay = 0;
	std::set<uint32_t> mKnownAreas;
	std::vector<uint32_t> mSoundSchemeIds; // At most 2
	bool mSoundGameOverSound = false;
	std::vector<uint32_t> mSoundGameOverStash; // 2 scheme ids
	bool mSoundGameInCombat = false;
	std::vector<uint32_t> mSoundGameCombatStash; // 2 scheme ids
	std::vector<QuestSave> mQuests;

	bool mInCombat = false;

	GameTime mRealTime;
	GameTime mGameTime;
	GameTime mAnimTime;
	std::vector<TimeEventSave> mTimeEvents;

	uint32_t mNextAnimId = 0;
	uint32_t mActiveGoalCount = 0;
	bool mAnimCatchup = false;
	uint32_t mAnimActionId = 0;
	std::vector<uint32_t> mUnkAnimVals;
	std::map<size_t, AnimSlotSave> mAnimSlots;

	static AnimSlotSave ReadAnimSlot(InputStream &in);
	static ObjRefSave ReadObjRef(InputStream &in);

	static void ReadSystemFooter(InputStream &in, const char *system);
};

max_stars_repo_name: edoipi/TemplePlus

max_issues_repo_name: edoipi/TemplePlus

max_forks_repo_name: edoipi/TemplePlus