text:
/* ode-initval/rk2simp.c
 *
 * Copyright (C) 2004 Tuomo Keskitalo
 *
 * 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.
 */

/* Runge-Kutta 2, Gaussian implicit. Also known as implicit midpoint rule.

   Non-linear equations solved by linearization, LU-decomposition
   and matrix inversion. For reference, see eg.

   Ascher, U.M., Petzold, L.R., Computer methods for ordinary
   differential and differential-algebraic equations, SIAM,
   Philadelphia, 1998.
 */

#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv.h>
#include <gsl/gsl_linalg.h>

#include "odeiv_util.h"

typedef struct
{
  double *Y1;
  double *y0;
  double *y0_orig;
  double *ytmp;
  double *dfdy;                 /* Jacobian */
  double *dfdt;                 /* time derivatives, not used */
  double *y_onestep;
  gsl_permutation *p;
}
rk2simp_state_t;

static void *
rk2simp_alloc (size_t dim)
{
  rk2simp_state_t *state =
    (rk2simp_state_t *) malloc (sizeof (rk2simp_state_t));

  if (state == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for rk2simp_state",
                      GSL_ENOMEM);
    }

  state->Y1 = (double *) malloc (dim * sizeof (double));

  if (state->Y1 == 0)
    {
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for Y1", GSL_ENOMEM);
    }

  state->y0 = (double *) malloc (dim * sizeof (double));

  if (state->y0 == 0)
    {
      free (state->Y1);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for y0", GSL_ENOMEM);
    }

  state->y0_orig = (double *) malloc (dim * sizeof (double));

  if (state->y0_orig == 0)
    {
      free (state->Y1);
      free (state->y0);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for y0_orig", GSL_ENOMEM);
    }

  state->ytmp = (double *) malloc (dim * sizeof (double));

  if (state->ytmp == 0)
    {
      free (state->Y1);
      free (state->y0);
      free (state->y0_orig);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for ytmp", GSL_ENOMEM);
    }

  state->dfdy = (double *) malloc (dim * dim * sizeof (double));

  if (state->dfdy == 0)
    {
      free (state->Y1);
      free (state->y0);
      free (state->y0_orig);
      free (state->ytmp);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for dfdy", GSL_ENOMEM);
    }

  state->dfdt = (double *) malloc (dim * sizeof (double));

  if (state->dfdt == 0)
    {
      free (state->Y1);
      free (state->y0);
      free (state->y0_orig);
      free (state->ytmp);
      free (state->dfdy);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for dfdt", GSL_ENOMEM);
    }

  state->y_onestep = (double *) malloc (dim * sizeof (double));

  if (state->y_onestep == 0)
    {
      free (state->Y1);
      free (state->y0);
      free (state->y0_orig);
      free (state->ytmp);
      free (state->dfdy);
      free (state->dfdt);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for y_onestep", GSL_ENOMEM);
    }

  state->p = gsl_permutation_alloc (dim);

  if (state->p == 0)
    {
      free (state->Y1);
      free (state->y0);
      free (state->y0_orig);
      free (state->ytmp);
      free (state->dfdy);
      free (state->dfdt);
      free (state->y_onestep);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for p", GSL_ENOMEM);
    }

  return state;
}

static int
rk2simp_step (double *y, rk2simp_state_t * state,
              const double h, const double t,
              const size_t dim, const gsl_odeiv_system * sys)
{
  /* Makes a Runge-Kutta 2nd order semi-implicit advance with step size h.
     y0 is initial values of variables y.

     The linearized semi-implicit equations to calculate are:

     Y1 = y0 + h/2 * (1 - h/2 * df/dy)^(-1) * f(t + h/2, y0)

     y = y0 + h * f(t + h/2, Y1)
   */

  const double *y0 = state->y0;
  double *Y1 = state->Y1;
  double *ytmp = state->ytmp;

  size_t i;
  int s, ps;

  gsl_matrix_view J = gsl_matrix_view_array (state->dfdy, dim, dim);

  /* First solve Y1.
     Calculate the inverse matrix (1 - h/2 * df/dy)^-1
   */

  /* Create matrix to J */

  s = GSL_ODEIV_JA_EVAL (sys, t, y0, state->dfdy, state->dfdt);

  if (s != GSL_SUCCESS)
    {
      return s;
    }

  gsl_matrix_scale (&J.matrix, -h / 2.0);
  gsl_matrix_add_diagonal(&J.matrix, 1.0);

  /* Invert it by LU-decomposition to invmat */

  s += gsl_linalg_LU_decomp (&J.matrix, state->p, &ps);

  if (s != GSL_SUCCESS)
    {
      return GSL_EFAILED;
    }

  /* Evaluate f(t + h/2, y0) */

  s = GSL_ODEIV_FN_EVAL (sys, t + 0.5 * h, y0, ytmp);

  if (s != GSL_SUCCESS)
    {
      return s;
    }

  /* Calculate Y1 = y0 + h/2 * ((1-h/2 * df/dy)^-1) ytmp */

  {
    gsl_vector_const_view y0_view = gsl_vector_const_view_array(y0, dim);
    gsl_vector_view ytmp_view = gsl_vector_view_array(ytmp, dim);
    gsl_vector_view Y1_view = gsl_vector_view_array(Y1, dim);

    s = gsl_linalg_LU_solve (&J.matrix, state->p,
                             &ytmp_view.vector, &Y1_view.vector);

    gsl_vector_scale (&Y1_view.vector, 0.5 * h);
    gsl_vector_add (&Y1_view.vector, &y0_view.vector);
  }

  /* And finally evaluation of f(t + h/2, Y1) and calculation of y */

  s = GSL_ODEIV_FN_EVAL (sys, t + 0.5 * h, Y1, ytmp);

  if (s != GSL_SUCCESS)
    {
      return s;
    }

  for (i = 0; i < dim; i++)
    {
      y[i] = y0[i] + h * ytmp[i];
    }

  return s;
}

static int
rk2simp_apply (void *vstate, size_t dim, double t, double h,
               double y[], double yerr[], const double dydt_in[],
               double dydt_out[], const gsl_odeiv_system * sys)
{
  rk2simp_state_t *state = (rk2simp_state_t *) vstate;

  size_t i;

  double *y0 = state->y0;
  double *y0_orig = state->y0_orig;
  double *y_onestep = state->y_onestep;

  /* Error estimation is done by step doubling procedure */

  DBL_MEMCPY (y0, y, dim);

  /* Save initial values in case of failure */
  DBL_MEMCPY (y0_orig, y, dim);

  /* First traverse h with one step (save to y_onestep) */
  DBL_MEMCPY (y_onestep, y, dim);

  {
    int s = rk2simp_step (y_onestep, state, h, t, dim, sys);

    if (s != GSL_SUCCESS)
      {
        return s;
      }
  }

  /* Then with two steps with half step length (save to y) */

  {
    int s = rk2simp_step (y, state, h / 2.0, t, dim, sys);

    if (s != GSL_SUCCESS)
      {
        /* Restore original y vector */
        DBL_MEMCPY (y, y0_orig, dim);
        return s;
      }
  }

  DBL_MEMCPY (y0, y, dim);

  {
    int s = rk2simp_step (y, state, h / 2.0, t + h / 2.0, dim, sys);

    if (s != GSL_SUCCESS)
      {
        /* Restore original y vector */
        DBL_MEMCPY (y, y0_orig, dim);
        return s;
      }
  }

  /* Derivatives at output */

  if (dydt_out != NULL)
    {
      int s = GSL_ODEIV_FN_EVAL (sys, t + h, y, dydt_out);

      if (s != GSL_SUCCESS)
        {
          /* Restore original y vector */
          DBL_MEMCPY (y, y0_orig, dim);
          return s;
        }
    }

  /* Error estimation */

  for (i = 0; i < dim; i++)
    {
      yerr[i] = 4.0 * (y[i] - y_onestep[i]) / 3.0;
    }

  return GSL_SUCCESS;
}

static int
rk2simp_reset (void *vstate, size_t dim)
{
  rk2simp_state_t *state = (rk2simp_state_t *) vstate;

  DBL_ZERO_MEMSET (state->Y1, dim);
  DBL_ZERO_MEMSET (state->y0, dim);
  DBL_ZERO_MEMSET (state->y0_orig, dim);
  DBL_ZERO_MEMSET (state->ytmp, dim);
  DBL_ZERO_MEMSET (state->dfdt, dim * dim);
  DBL_ZERO_MEMSET (state->dfdt, dim);
  DBL_ZERO_MEMSET (state->y_onestep, dim);

  return GSL_SUCCESS;
}

static unsigned int
rk2simp_order (void *vstate)
{
  rk2simp_state_t *state = (rk2simp_state_t *) vstate;
  state = 0;                    /* prevent warnings about unused parameters */
  return 2;
}

static void
rk2simp_free (void *vstate)
{
  rk2simp_state_t *state = (rk2simp_state_t *) vstate;
  free (state->Y1);
  free (state->y0);
  free (state->y0_orig);
  free (state->ytmp);
  free (state->dfdy);
  free (state->dfdt);
  free (state->y_onestep);
  gsl_permutation_free (state->p);
  free (state);
}

static const gsl_odeiv_step_type rk2simp_type = {
  "rk2simp",                    /* name */
  0,                            /* can use dydt_in? */
  1,                            /* gives exact dydt_out? */
  &rk2simp_alloc,
  &rk2simp_apply,
  &rk2simp_reset,
  &rk2simp_order,
  &rk2simp_free
};

const gsl_odeiv_step_type *gsl_odeiv_step_rk2simp = &rk2simp_type;


text:
#pragma once

#include <winrt\Windows.Foundation.h>
#include <d3d11.h>
#include <DirectXMath.h>
#include <gsl\gsl>
#include "DrawableGameComponent.h"
#include "MatrixHelper.h"
#include "VectorHelper.h"
#include "BasicMaterial.h"

namespace Library
{
	class Mesh;

	class ProxyModel final : public DrawableGameComponent
	{
		RTTI_DECLARATIONS(ProxyModel, DrawableGameComponent)

	public:
		ProxyModel(Game& game, const std::shared_ptr<Camera>& camera, const std::string& modelFileName, float scale = 1.0f);
		ProxyModel(const ProxyModel&) = delete;
		ProxyModel(ProxyModel&&) = default;
		ProxyModel& operator=(const ProxyModel&) = delete;
		ProxyModel& operator=(ProxyModel&&) = default;
		~ProxyModel() = default;

		const DirectX::XMFLOAT3& Position() const;
		const DirectX::XMFLOAT3& Direction() const;
		const DirectX::XMFLOAT3& Up() const;
		const DirectX::XMFLOAT3& Right() const;

		DirectX::XMVECTOR PositionVector() const;
		DirectX::XMVECTOR DirectionVector() const;
		DirectX::XMVECTOR UpVector() const;
		DirectX::XMVECTOR RightVector() const;

		bool& DisplayWireframe();

		void SetPosition(float x, float y, float z);
		void SetPosition(DirectX::FXMVECTOR position);
		void SetPosition(const DirectX::XMFLOAT3& position);

		void ApplyRotation(DirectX::CXMMATRIX transform);
		void ApplyRotation(const DirectX::XMFLOAT4X4& transform);

		void SetColor(const DirectX::XMFLOAT4& color);

		virtual void Initialize() override;
		virtual void Update(const GameTime& gameTime) override;
		virtual void Draw(const GameTime& gameTime) override;

	private:
		DirectX::XMFLOAT4X4 mWorldMatrix{ MatrixHelper::Identity };
		DirectX::XMFLOAT3 mPosition{ Vector3Helper::Zero };
		DirectX::XMFLOAT3 mDirection{ Vector3Helper::Forward };
		DirectX::XMFLOAT3 mUp{ Vector3Helper::Up };
		DirectX::XMFLOAT3 mRight{ Vector3Helper::Right };
		BasicMaterial mMaterial;
		std::string mModelFileName;
		float mScale;
		winrt::com_ptr<ID3D11Buffer> mVertexBuffer;
		winrt::com_ptr<ID3D11Buffer> mIndexBuffer;
		std::uint32_t mIndexCount{ 0 };
		bool mDisplayWireframe{ true };
		bool mUpdateWorldMatrix{ true };
		bool mUpdateMaterial{ true };
	};
}

ext: h


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

/*
 * rng.h
 *
 *  Created on: Jul 3, 2013
 *      Author: marc
 */

#ifndef RNG_H_
#define RNG_H_

#include <climits> // ULONG_MAX
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <iostream>

// Seeds for GNU GSL random number generators
typedef unsigned long int gslseed_t;

class RandomNumberGenerator
{
public:
    /*
    RandomNumberGenerator()
    : mRngType( gsl_rng_default ), mSeed( time(0) ) {
    	initializeRNG();
    }*/

    RandomNumberGenerator(const gslseed_t& seed)
        : mRngType( gsl_rng_default ), mSeed( seed )
    {
        initializeRNG();
    }

    RandomNumberGenerator(const RandomNumberGenerator& rhs)
    {
        mRngType = rhs.mRngType;
        mRng = gsl_rng_clone( rhs.mRng );
        mSeed = rhs.randomSeed();
        gsl_rng_set (mRng, mSeed);
    }

    ~RandomNumberGenerator()
    {
        gsl_rng_free( mRng );
    }

    // Random uniform number between 0 and 1.
    double uniform() const
    {
        return gsl_rng_uniform( mRng );
    }

    // Random uniform number between min and max
    double uniform(const double min, const double max) const
    {
        const double u = uniform();
        return ( min + u * (max - min) );
    }

    double gaussian(const double sigma) const
    {
        return ( gsl_ran_gaussian (mRng, sigma) );
    }

    // Beta distribution. Random x between 0 and 1.
    double beta(const double alpha, const double beta) const
    {
        return ( gsl_ran_beta(mRng, alpha, beta) );
    }

    gslseed_t randomSeed() const
    {
    	// I wrote the following after looking into the GSL code itself.
        const unsigned long int offset = mRng->type->min;
        const unsigned long int range = mRng->type->max - offset;
        gslseed_t rseed = 1 + gsl_rng_uniform_int(mRng, range-2);
        //std::cerr << "generated random seed = " << rseed << std::endl;
        return rseed;
    }

protected:
    void initializeRNG()
    {
        mRng = gsl_rng_alloc (mRngType);
        gsl_rng_set (mRng, mSeed);
    }

private:
    // Prevent copying and assignment
    //RandomNumberGenerator(const RandomNumberGenerator&);
    void operator=(const RandomNumberGenerator&);

    const gsl_rng_type* mRngType;
    gslseed_t			mSeed;
    gsl_rng* 			mRng;
};

#endif /* RNG_H_ */

max_stars_repo_name: marcnormandin/ParticleSwarmOptimization

max_issues_repo_name: marcnormandin/ParticleSwarmOptimization

max_forks_repo_name: marcnormandin/ParticleSwarmOptimization


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

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_movstat.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_statistics.h>

typedef double funcacc_type_t;
typedef funcacc_type_t ringbuf_type_t;

#include "ringbuf.c"

typedef struct
{
  funcacc_type_t *window; /* linear array for current window */
  ringbuf *rbuf;          /* ring buffer storing current window */
} funcacc_state_t;

static size_t
funcacc_size(const size_t n)
{
  size_t size = 0;

  size += sizeof(funcacc_state_t);
  size += n * sizeof(funcacc_type_t);
  size += ringbuf_size(n);

  return size;
}

static int
funcacc_init(const size_t n, void * vstate)
{
  funcacc_state_t * state = (funcacc_state_t *) vstate;

  state->window = (funcacc_type_t *) ((unsigned char *) vstate + sizeof(funcacc_state_t));
  state->rbuf = (ringbuf *) ((unsigned char *) state->window + n * sizeof(funcacc_type_t));

  ringbuf_init(n, state->rbuf);

  return GSL_SUCCESS;
}

static int
funcacc_insert(const funcacc_type_t x, void * vstate)
{
  funcacc_state_t * state = (funcacc_state_t *) vstate;

  /* add new element to ring buffer */
  ringbuf_insert(x, state->rbuf);

  return GSL_SUCCESS;
}

static int
funcacc_delete(void * vstate)
{
  funcacc_state_t * state = (funcacc_state_t *) vstate;

  if (!ringbuf_is_empty(state->rbuf))
    ringbuf_pop_back(state->rbuf);

  return GSL_SUCCESS;
}

static int
funcacc_get(void * params, funcacc_type_t * result, const void * vstate)
{
  const funcacc_state_t * state = (const funcacc_state_t *) vstate;
  gsl_movstat_function *f = (gsl_movstat_function *) params;
  size_t n = ringbuf_copy(state->window, state->rbuf);

  *result = GSL_MOVSTAT_FN_EVAL(f, n, state->window);

  return GSL_SUCCESS;
}

static const gsl_movstat_accum func_accum_type =
{
  funcacc_size,
  funcacc_init,
  funcacc_insert,
  funcacc_delete,
  funcacc_get
};

const gsl_movstat_accum *gsl_movstat_accum_userfunc = &func_accum_type;

ext: c

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
#ifndef BIO_RUN_MATCH_H_
#define BIO_RUN_MATCH_H_

#include "bio/defs.h"
#include "bio/common.h"
#include "bio/match_hit.h"
#include "bio/sequence.h"
#include "bio/biobase_filter.h"
#include "bio/pssm_match.h"
#include "bio/phylogenetics.h"
#include "bio/biobase_db.h"
#include "bio/biobase_data_traits.h"

#include <boost/filesystem/path.hpp>

#include <vector>
#include <set>

#include <gsl/gsl_sf_log.h>
#include <gsl/gsl_sf_exp.h>

namespace boost { namespace program_options {
    //forward decl
    class options_description;
} }

BIO_NS_START

enum ScoreAlgorithm
{
    OTT_SCORE_ALGORITHM,
    BAYESIAN_SCORE_ALGORITHM
};

struct MatchParams
{
    MatchParams(
        float_t t,
        ScoreAlgorithm alg = OTT_SCORE_ALGORITHM,
        bool or_better = false)
        : threshold( t )
        , score_algorithm( alg )
        , or_better( or_better )
    {
    }

    float_t threshold;
    ScoreAlgorithm score_algorithm;
    bool or_better;
};

struct MatchResults
{
    MatchResults();

    MatchResults(
        TableLink link,
        Hit result);

    TableLink link;
    Hit result;
    size_t number;

    bool operator<(const MatchResults & rhs) const
    {
        if (result == rhs.result)
        {
            if (link == rhs.link)
            {
                return number < rhs.number;
            }
            else
            {
                return link < rhs.link;
            }
        }
        else
        {
            return result < rhs.result;
        }
    }

protected:
    friend class boost::serialization::access;
    template<class Archive>
    void serialize(Archive & ar, const unsigned int version)
    {
        ar & link;
        ar & result;
        ar & number;
    }
};
typedef std::vector<MatchResults> match_result_vec_t;
typedef boost::shared_ptr< match_result_vec_t > match_result_vec_ptr_t;
void sort_by_position(match_result_vec_t & matches);
std::ostream & operator<<(std::ostream & os, const MatchResults & match);

/** Compare based initially on position, then link. */
struct MatchResultsPositionLessThan
{
    bool
    operator()(const MatchResults & h1, const MatchResults & h2) const
    {
        return
            h1.result.position < h2.result.position
            ||
            (
                h1.result.position == h2.result.position
                &&
                h1.link < h2.link
            );
    }
};

inline
std::ostream &
operator<<(std::ostream & os, const MatchResults & details) {
    os
        << details.link << ","
        << details.result;
    return os;
}

/**
Score all the pssm on the sequence between match_seq_begin and match_seq_end and put the results in result_insert_it.
Returns an estimate that the pssm binds the sequence.
*/
template <class Pssm, class ResultInsIt>
float_t
score_pssm(
    const Pssm & pssm,
    seq_t::const_iterator match_seq_begin,
    seq_t::const_iterator match_seq_end,
    const MatchParams & params,
    ResultInsIt & result_insert_it)
{
    Scorers scorers(match_seq_begin, match_seq_end);

    //const hit_vec_t & hit_results = scorers.bayesian_scores.get_result(pssm);
    const hit_vec_t & hit_results =
        OTT_SCORE_ALGORITHM == params.score_algorithm
            ? scorers.ott_normalised_scores.get_result(pssm)
            : ( params.or_better
                ? scorers.bayesian_scores.get_result(pssm)
                : scorers.bayesian_or_better_scores.get_result(pssm));

    float_t p_does_not_bind = 1.0;
    for (hit_vec_t::const_iterator i = hit_results.begin();
        hit_results.end() != i;
        ++i)
    {
        if (i->score > params.threshold)
        {
            p_does_not_bind *= ( float_t( 1.0 ) - i->score );

            MatchResults result(pssm.get_link(), *i);
            assert(result.result.position < (int) (match_seq_end - match_seq_begin)); //make sure the position is not too high
            *result_insert_it++ = result;
        }
    }
    return float_t( 1.0 ) - p_does_not_bind;
}

/** Score all the pssms between begin and end and put the results in result_insert_it. */
template <class PssmIt, class ResultInsIt>
size_t
score_pssms(
    PssmIt pssm_begin,
    PssmIt pssm_end,
    seq_t::const_iterator match_seq_begin,
    seq_t::const_iterator match_seq_end,
    const MatchParams & params,
    ResultInsIt & result_insert_it)
{
    size_t num_pssms_matched = 0;
    for ( ; pssm_begin != pssm_end; ++pssm_begin)
    {
        try
        {
            score_pssm(
                *(pssm_begin->second.get()),
                match_seq_begin,
                match_seq_end,
                params,
                result_insert_it);

            ++num_pssms_matched;
        }
        catch (const std::exception & ex)
        {
            std::cerr << "Could not score " << pssm_begin->second->get_name() << ": " << ex.what() << std::endl;
        }
        catch (const std::string & ex)
        {
            std::cerr << "Could not score " << pssm_begin->second->get_name() << ": " << ex << std::endl;
        }
        catch (const char * ex)
        {
            std::cerr << "Could not score " << pssm_begin->second->get_name() << ": " << ex << std::endl;
        }
        catch (...)
        {
            std::cerr << "Could not score " << pssm_begin->second->get_name() << std::endl;
        }
    }
    return num_pssms_matched;
}

template <class SeqIt, class ResultInsIt, class ConsFilter, class MatrixFilter>
void
pssm_match(
    SeqIt seq_begin,
    SeqIt seq_end,
    float_t threshold,
    ScoreAlgorithm algorithm,
    bool use_or_better,
    ConsFilter consensus_filter,
    MatrixFilter matrix_filter,
    ResultInsIt result_inserter)
{
    typedef boost::filter_iterator<MatrixFilter, Matrix::map_t::const_iterator> matrix_filter_it;
    typedef boost::filter_iterator<ConsFilter, Site::map_t::const_iterator> site_filter_it;

    const matrix_filter_it matrices_begin(matrix_filter, BiobaseDb::singleton().get_matrices().begin(), BiobaseDb::singleton().get_matrices().end());
    const matrix_filter_it matrices_end(matrix_filter, BiobaseDb::singleton().get_matrices().end(), BiobaseDb::singleton().get_matrices().end());

    const site_filter_it sites_begin(consensus_filter, BiobaseDb::singleton().get_sites().begin(), BiobaseDb::singleton().get_sites().end());
    const site_filter_it sites_end(consensus_filter, BiobaseDb::singleton().get_sites().end(), BiobaseDb::singleton().get_sites().end());

    if (matrices_begin == matrices_end && sites_begin == sites_end)
    {
        throw std::logic_error( "No PSSMs match filter" );
    }

    MatchParams params(threshold, algorithm, use_or_better);
    {
        //std::cout << "Scoring matrices" << std::endl;
        //const size_t num_matched =
            score_pssms(
                matrices_begin,
                matrices_end,
                seq_begin,
                seq_end,
                params,
                result_inserter);
//                ostream_iterator<MatchResults>(cout, "\n"));
        //std::cout << "Scored " << num_matched << " matrices" << std::endl;
    }

    {
        //std::cout << "Scoring sites" << std::endl;
        //const size_t num_matched =
            score_pssms(
                sites_begin,
                sites_end,
                seq_begin,
                seq_end,
                params,
                result_inserter);
//                ostream_iterator<MatchResults>(cout, "\n"));
        //std::cout << "Scored " << num_matched << " sites" << std::endl;
    }
}

void
pssm_match(
    const seq_t & match_seq,
    float_t threshold,
    ScoreAlgorithm algorithm,
    const boost::filesystem::path & file,
    const std::string & title,
    bool show_labels);

template <typename HitIt>
float_t
estimate_binding_prob(
    HitIt hit_begin,
    HitIt hit_end)
{
    float_t product = 1.0;
    for (HitIt hit = hit_begin; hit_end != hit; ++hit)
    {
        product *= float_t( 1.0f - hit->result.score );
    }
    return float_t(1.0 - product);
}

template <typename Exponent>
double
power(double x, Exponent y)
{
    return
        0 == x
            ? 0
            : gsl_sf_exp(y * gsl_sf_log(x));
}

/** Adjusts hits for occurence in phylogenetically conserved sequence. */
template <typename ResultIt>
void
adjust_hits(
    ResultIt results_begin,
    ResultIt results_end,
    seq_t::const_iterator seq_begin,
    seq_t::const_iterator seq_end,
    float_t threshold,
    ScoreAlgorithm algorithm)
{
    const MatchParams params(threshold, algorithm);

    //maintain a set of already adjusted PSSMs
    std::set<TableLink> already_adjusted;

    //for each result
    for (ResultIt result = results_begin; results_end != result; ++result)
    {
        //check we haven't checked this PSSM before
        if (already_adjusted.end() != already_adjusted.find(result->link))
        {
            //we have so ignore it
            continue;
        }

        //we need to estimate likelihood that it binds in phylogenetic sequence so run PSSM over sequence.
        match_result_vec_t phylo_hits;
        std::insert_iterator<bio::match_result_vec_t> inserter(phylo_hits, phylo_hits.begin());
        switch (result->link.table_id)
        {
        case MATRIX_DATA:
            score_pssm(
                *BiobaseDb::singleton().get_entry<MATRIX_DATA>(result->link),
                seq_begin,
                seq_end,
                params,
                inserter);
            break;

        case SITE_DATA:
            score_pssm(
                *BiobaseDb::singleton().get_entry<SITE_DATA>(result->link),
                seq_begin,
                seq_end,
                params,
                inserter);
            break;

        default:
            throw std::logic_error( "Unknown pssm type" );
        }

        const float_t binding_prob = estimate_binding_prob(phylo_hits.begin(), phylo_hits.end());

        //adjust all the hits for this PSSM
        for (ResultIt to_adjust = result; results_end != to_adjust; ++to_adjust)
        {
            if (to_adjust->link == result->link)
            {
                //to_adjust->result.score *= power(binding_prob, 5.0 * phylo_seq.conservation * phylo_seq.conservation);
                to_adjust->result.score *= binding_prob;
            }
        }

        //add to the list of already adjusted PSSMs
        already_adjusted.insert(result->link);
    }
}

/** Adjusts hits for occurence in phylogenetically conserved sequence. */
template <
    typename ResultIt,
    typename SeqRange >
void
adjust_hits_for_phylo_sequences(
    ResultIt results_begin,
    ResultIt results_end,
    const SeqRange & sequences,
    float_t threshold,
    ScoreAlgorithm algorithm)
{
    unsigned num_seqs = 0;
    BOOST_FOREACH( const typename boost::range_value< SeqRange >::type & seq, sequences )
    {
        adjust_hits(
            results_begin,
            results_end,
            seq.begin(),
            seq.end(),
            threshold,
            algorithm );
        ++num_seqs;
    }
    std::cout << "Adjusted hits for " << num_seqs << " sequences" << std::endl;

    for( ; results_begin != results_end; ++results_begin )
    {
        //adjust hits - raise to 1/#seqs
        results_begin->result.score =
            float_t(
                ( 0.0 == results_begin->result.score ) ? 0.0 : exp( log( results_begin->result.score ) / ( num_seqs + 1 ) ) );
    }
}

/** Adjust hits by power. */
struct RaiseHitToPower
{
    double power;

    RaiseHitToPower(double power);

    void operator()(MatchResults & hit) const;
};

/** Is a hit above a threshold? */
struct HitAboveThreshold
{
    float_t threshold;

    HitAboveThreshold(float_t threshold);

    bool operator()(const MatchResults & results) const;
};

BIO_NS_END

#endif //BIO_RUN_MATCH_H_

max_stars_repo_name: JohnReid/biopsy

max_issues_repo_name: JohnReid/biopsy

max_forks_repo_name: JohnReid/biopsy


text:
/* ============================================================ *
 * lensing_3rd.h						*
 * Martin Kilbinger, Liping Fu 2010.				*
 * ============================================================ */

#ifndef __LENSING_3RD_H
#define __LENSING_3RD_H

#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <fftw3.h>

#include <gsl/gsl_integration.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_plain.h>
#include <gsl/gsl_monte_vegas.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf_bessel.h>

#include "cosmo.h"
#include "nofz.h"
#include "lensing.h"
#include "errorlist.h"
#include "maths_base.h"

#define tenoverseven   1.428571428571
#define fouroverseven  0.571428571429

/* Minimum scale factor where a non-linear scale can be defined. *
 * See {a,b,c}scocou.						 */
#define A_NL_MIN 0.02

#define s2_min     0.1
#define s2_max     1.0e6
#define N_s2       50
#define epsilon0   1.0e-2
#define N_EFF_MIN -2.0

/* more ce_xyz in smith2.h */
#define lensing_3rd_base         -2400
#define lensing_3rd_wrongmode    -1 + lensing_3rd_base
#define lensing_3rd_rootbracket  -4 + lensing_3rd_base
#define lensing_3rd_slc          -5 + lensing_3rd_base

typedef enum {fgauss=0, fpoly=1, ftophat=2, fdelta=3, fxip=4, fxim=5, fall=6} filter_t;

typedef enum {PT=0, SCOCOU=1, GM12} bispmode_t;
#define sbispmode_t(i) (   \
  i==PT     ? "PT" :       \
  i==SCOCOU ? "scocou01" : \
  i==GM12   ? "GM12" :     \
  "")
#define Nbispmode_t 3

/* Intrinsic 3rd-order alignment model */
typedef enum {ia_3rd_none, ia_3rd_S08} ia_3rd_t;
#define sia_3rd_t(i) ( \
  i==ia_3rd_none ? "none" : \
  i==ia_3rd_S08  ? "S08" : \
  "")
#define Nia_3rd_t 2

/* Bit-coded IA terms */
typedef enum {ia_3rd_undef, ia_GGI_GII_III, ia_only_GGI, ia_only_GII, ia_only_III} ia_3rd_terms_t;
#define sia_3rd_terms_t(i) (      \
  i==ia_3rd_undef     ? "undef" :   \
  i==ia_GGI_GII_III   ? "GGI_GII_III" :   \
  i==ia_only_GGI      ? "only_GGI" : \
  i==ia_only_GII      ? "only_GII" : \
  i==ia_only_III      ? "only_III" : \
  "")
#define Nia_3rd_terms_t 5

/* Source-lens clustering */
typedef enum {slc_none, slc_FK13} slc_t;
#define sslc_t(i) ( \
  i==slc_none      ? "none" : \
  i==slc_FK13      ? "slc_FK13" : \
  "")
#define Nslc_t 2

typedef enum {kkk=0, kkg=1, kgg=2, ggg=3} bispfield_t;

typedef struct {

  /* Lensing, including basic cosmology and redshift distribution)(s) */
  cosmo_lens *lens;

  /* Intrinsic alignment parameters */
  ia_3rd_t ia;
  ia_3rd_terms_t ia_terms;
  double A_GGI, theta_GGI, A_GII, theta_GII;

  /* Source-lens clustering */
  slc_t slc;
  double b_slc, gamma_slc;

  /* ============================================================ *
   * Precomputed stuff (at the moment only one redshift-bin).     *
   * ============================================================ */

  interTable *k_NL, *n_eff;
  double scale_NL_amin;
  interTable2D **B_kappa[3];

  bispmode_t bispmode;

} cosmo_3rd;

typedef struct {
  double r1, r2;
  cosmo_3rd *self;
} cosmo3ANDdoubleANDdouble;

typedef struct {
  double r1, r2;
  cosmo_3rd *self;
  int i, n_bin[3];
} cosmo3ANDtomo;

typedef struct {
  cosmo_3rd *self;
  double R[3];
  filter_t wfilter;
  int n_bin[3];
  int m;
  error **err;
} cosmo3ANDmorestuff;

typedef struct {
   cosmo_3rd *self;
   error **err;
   double a1, a2, f1, f2, R;
} cosmo3SLC;

cosmo_3rd *init_parameters_3rd(double OMEGAM, double OMEGAV, double W0_DE, double W1_DE,
			       double *W_POLY_DE, int N_POLY_DE,
			       double H100, double OMEGAB, double OMEGANUMASS,
			       double NEFFNUMASS, double NORM, double NSPEC,
			       int Nzbin, const int *Nnz, const nofz_t *nofz, double *par_nz,
			       nonlinear_t NONLINEAR, transfer_t TRANSFER,
			       growth_t GROWTH, de_param_t DEPARAM,
			       norm_t normmode,
			       ia_t IA, ia_terms_t IA_TERMS, double A_IA,
			       bispmode_t BISPMODE,
			       ia_3rd_t IA_3RD, ia_3rd_terms_t IA_3RD_TERMS, double A_GGI, double theta_GGI,
			       double A_GII, double theta_GII,
			       slc_t slc, double b_slc, double gamma_slc,
			       error **err);

void consistency_parameters_3rd(const cosmo_3rd *self, error **err);
cosmo_3rd *copy_parameters_3rd_only(cosmo_3rd *source, error **err);
void updateFrom_3rd(cosmo_3rd *avant, cosmo_3rd *apres, error **err);
cosmo_3rd *set_cosmological_parameters_to_default_lens_3rd(error **err);
void read_cosmological_parameters_lens_3rd(cosmo_3rd **self, FILE *F, error **err);
void free_parameters_3rd(cosmo_3rd **self);
void dump_param_3rd(cosmo_3rd* self, FILE *F, error **err);

double dcub(double a);
double n_eff_one(cosmo *self, double k, error **err);
double n_eff(cosmo_3rd *, double, int, error **);
double Q3(double, error **err);
double temp_NL(double, double, cosmo *, error **);
double scale_NL(cosmo_3rd *, double, error **);

double ascocou(cosmo_3rd *, double, double, error **);
double bscocou(cosmo_3rd *, double, double, error **);
double cscocou(cosmo_3rd *, double, double, error **);

double int_for_B_kappa_bar0(double a, void *intpar, error **err);
double int_for_B_kappa_bar1(double a, void *intpar, error **err);
double int_for_B_kappa_bar2(double a, void *intpar, error **err);
double int_for_B_kappa_bar3(double a, void *intpar, error **err);

double F2eff(cosmo_3rd *, double a, double k1, double k2, double cosphi, error **);
double F2bar(int, double, double, error **);
double F2(int, double, double, error **);
double F2cos(int, double, error **err);

double Q_123(cosmo_3rd *, double, double, double, double, error **);
double bb(cosmo_3rd *, double, double, double, int i_bin, int j_bin, int k_bin, error **);
double B_kappa_bar(cosmo_3rd *self, double s1, double s2, int abc, int i_bin, int j_bin, int k_bin, error **err);

double hept_rtbis(double (*func)(double,double,cosmo*,error**), double, double,
		  double, double, cosmo *, error **);
double B_delta(cosmo_3rd *self, double k1, double k2, double cosbeta, double a,
	       error **err);

double Uhat_one(double x, filter_t wfilter);
double Uhat(double eta, filter_t wfilter);
void permute3(double *x, int offset);
double int_for_map3_3d(double x[], size_t dim, void *intpar);
double map3_perm(cosmo_3rd *self, double R[3], int i_bin, int j_bin, int k_bin, filter_t wfilter, error **err);
double map3(cosmo_3rd *self, double R[3], int i_bin, int j_bin, int k_bin, filter_t wfilter, error **err);

double int_for_E_GGI(cosmo_lens *self, double zs1, double zs2, double zl, error **err);
double E_GGI(cosmo_lens *self, error **err);
double map3_GGI(cosmo_3rd *self, double theta, error **err);

double int_for_E_GII(cosmo_lens *self, double zs, double zl, error **err);
double E_GII(cosmo_lens *self, error **err);
double map3_GII(cosmo_3rd *self, double theta, error **err);

double map3_SLC_t1(cosmo_3rd *self, double R, int n_bin, error **err);
double bias_SLC(cosmo_3rd *self, double a, error **err);
double int_for_Q_mc(double x[], size_t dim, void *intpar);
double Q_mc(cosmo_3rd *self, double a1, double a2, double f1, double f2, double R, double *abserr,
	    double *chisqr, error **err);

double lensing_signal_3rd(cosmo_3rd *self, double theta[3], int i_bin, int j_bin, int k_bin, error **err);
void fill_dmm_map3gauss_diag(cosmo_3rd *self, double *data_minus_model, int start, const double *data,
			     int Nzbin, int Ntheta, double *theta, error **err);
void fill_dmm_map3gauss(cosmo_3rd *self, double *data_minus_model, int start, const double *data,
			int Ntheta, double *theta, error **err);
double chi2_lensing_3rd(cosmo_3rd *self, datcov *dc, const cosebi_info_t *cosebi_info, error **err);

#endif

max_stars_repo_name: danielgruen/ccv

max_issues_repo_name: danielgruen/ccv

max_forks_repo_name: danielgruen/ccv


text:
#pragma once

#include <gsl/gsl>

/*
	Helper class for reading structured binary data from a
	memory stream.
*/
class BinaryReader {
public:

	explicit BinaryReader(gsl::span<uint8_t> data) : mData(data) {}

	template<typename T>
	T Read() {
		Expects(mData.size() >= sizeof(T));
		auto result{ *reinterpret_cast<T*>(&mData[0]) };
		mData = mData.subspan(sizeof(T));
		return result;
	}

	std::string ReadFixedString(size_t length) {
		Expects(mData.size() >= (int) length);
		std::string result(length, '\0');
		memcpy(&result[0], &mData[0], length);
		result.resize(result.length());
		mData = mData.subspan(length);
		return result;
	}

	bool AtEnd() const {
		return mData.size() == 0;
	}

private:
	gsl::span<uint8_t> mData;
};

max_stars_repo_name: edoipi/TemplePlus

max_issues_repo_name: edoipi/TemplePlus

max_forks_repo_name: edoipi/TemplePlus


text:
/*
** Implementation of LISA algorithm
** for statistical inference of fMRI images
**
** 1st level analysis using GLM (general linear model) with pre-coloring
**
** G.Lohmann, April 2017
*/
#include "viaio/Vlib.h"
#include "viaio/file.h"
#include "viaio/mu.h"
#include "viaio/option.h"
#include "viaio/os.h"
#include <viaio/VImage.h>

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_permutation.h>

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#ifdef _OPENMP
#include <omp.h>
#endif /*_OPENMP*/

typedef struct TrialStruct {
  int   id;
  float onset;
  float duration;
  float height;
} Trial;

extern void VIsolatedVoxels(VImage src,float threshold);
extern void VHistogram(gsl_histogram *histogram,VString filename);
extern void VCheckImage(VImage src);
extern void FDR(VImage src,VImage dest,gsl_histogram *nullhist,gsl_histogram *realhist,double);
extern double ttest1(double *data1,int n);
extern void ImageStats(VImage src,double *,double *,double *hmin,double *hmax);
extern Trial *ReadDesign(VStringConst designfile,int *numtrials,int *nevents);
extern gsl_matrix *VCreateDesign(int ntimesteps,int nevents,int,VBoolean,gsl_matrix *);
extern void VHemoModel(Trial *trial,int ntrials,int nevents,int ntimesteps,double tr,int,VBoolean,gsl_matrix *,gsl_matrix *);
extern Trial *CopyTrials(Trial *trial,int numtrials);
extern void VGLM(gsl_matrix *Data,gsl_matrix *X,gsl_matrix *XInv,gsl_vector *con,VImage map,VImage zmap);
extern Trial *ConcatenateTrials(Trial **trial,int *numtrials,float *run_duration,int dlists,int sumtrials);

extern double VImageVar(VImage src);
extern void VImageCount(VImage src);
extern void VBilateralFilter(VImage src,VImage,int radius,double var1,double var2,int);
extern void VGetHistRange(VImage src,double *hmin,double *hmax);
extern void VZScale(VImage src,float,float stddev);
extern float VGetMode(VImage src);
extern void GlobalMean(gsl_matrix *Data,gsl_matrix *covariates,int column);
extern gsl_matrix *VReadCovariates(VString cfile,VBoolean normalize);

extern VImage VoxelMap(VAttrList list);
extern gsl_matrix *VReadImageData(VAttrList *list,int nlists);
extern void VGetTimeInfos(VAttrList *list,int nlists,double *mtr,float *run_duration);
extern void VRowNormalize(gsl_matrix *Data);
extern void CheckTrialLabels(Trial *trial,int numtrials);
extern void HistoUpdate(VImage,gsl_histogram *);
extern void PlotDesign(gsl_matrix *X,double tr,VString filename);

void XCheckImage(VImage src,char *filename)
{
  VAttrList out_list = VCreateAttrList();
  VAppendAttr(out_list,"image",NULL,VImageRepn,src);
  FILE *out_file = fopen(filename,"w");
  VWriteFile (out_file, out_list);
}

/* shuffle each run separately to ensure exchangebility, concatenate individual permtables */
int **genperm(gsl_rng *rx,int *numtrials,int sumtrials,int dlists,int numperm)
{
  int i,j,k;

  int **permtable = (int **) VCalloc(numperm,sizeof(int *));

  gsl_permutation **perm = (gsl_permutation **) VCalloc(dlists,sizeof(gsl_permutation *));
  for (k=0; k<dlists; k++) {
    perm[k] = gsl_permutation_alloc((size_t)numtrials[k]);
    gsl_permutation_init (perm[k]);
  }

  for (i = 0; i < numperm; i++) {
    permtable[i] = (int *) VCalloc(sumtrials,sizeof(int));
    int jj=0;
    for (k=0; k<dlists; k++) {
      gsl_ran_shuffle (rx, perm[k]->data,numtrials[k],sizeof(size_t));
      for (j=0; j<numtrials[k]; j++) {
	permtable[i][j+jj] = perm[k]->data[j] + jj;
      }
      jj += numtrials[k];
    }
  }

  for (k=0; k<dlists; k++) {
    gsl_permutation_free(perm[k]);
  }
  return permtable;
}

VDictEntry HemoDict[] = {
  { "gamma_0", 0 },
  { "gamma_1", 1 },
  { "gamma_2", 2 },
  { "gauss", 3 },
  { NULL }
};

int main (int argc, char *argv[])
{
  static VArgVector in_files;
  static VArgVector des_files;
  static VString  cova_filename="";
  static VString  out_filename="";
  static VString  plot_filename="";
  static VString  mask_filename="";
  static VShort   hemomodel = 0;
  static VBoolean firstcol = TRUE;
  static VArgVector contrast;
  static VFloat   alpha = 0.05;
  static VShort   radius = 2;
  static VFloat   rvar = 2.0;
  static VFloat   svar = 2.0;
  static VShort   numiter = 2;
  static VBoolean cleanup = TRUE;
  static VBoolean demean = TRUE;
  static VBoolean verbose = FALSE;
  static VBoolean globalmean = FALSE;
  static VShort   numperm = 5000;
  static VLong    seed = 99402622;
  static VShort   nproc = 0;
  static VOptionDescRec options[] = {
    {"in", VStringRepn, 0, & in_files, VRequiredOpt, NULL,"Input files" },
    {"out", VStringRepn, 1, & out_filename, VRequiredOpt, NULL,"Output file" },
    {"design", VStringRepn, 0, & des_files, VRequiredOpt, NULL,"Design files (1st level)" },
    {"contrast", VFloatRepn, 0, (VPointer) &contrast, VRequiredOpt, NULL, "Contrast vector"},
    {"nuisance", VStringRepn,  1, & cova_filename, VOptionalOpt, NULL,"Nuisance regressors" },
    {"demean",VBooleanRepn,1,(VPointer) &demean,VOptionalOpt,NULL,"Whether to subtract mean in nuisance regressors"},
    {"hemo", VShortRepn, 1, (VPointer) &hemomodel, VOptionalOpt, HemoDict,"Hemodynamic model" },
    {"col1", VBooleanRepn, 1, (VPointer) &firstcol, VOptionalOpt, NULL,"Whether to add a constant first column" },
    {"alpha",VFloatRepn,1,(VPointer) &alpha,VOptionalOpt,NULL,"FDR significance level"},
    {"perm",VShortRepn,1,(VPointer) &numperm,VOptionalOpt,NULL,"Number of permutations"},
    {"seed",VLongRepn,1,(VPointer) &seed,VOptionalOpt,NULL,"Seed for random number generation"},
    {"plotdesign", VStringRepn, 1, & plot_filename, VOptionalOpt, NULL,"Filename for plotting design matrix X" },
    {"radius",VShortRepn,1,(VPointer) &radius,VOptionalOpt,NULL,"Bilateral parameter (radius in voxels)"},
    {"rvar",VFloatRepn,1,(VPointer) &rvar,VOptionalOpt,NULL,"Bilateral parameter (radiometric)"},
    {"svar",VFloatRepn,1,(VPointer) &svar,VOptionalOpt,NULL,"Bilateral parameter (spatial)"},
    {"filteriterations",VShortRepn,1,(VPointer) &numiter,VOptionalOpt,NULL,"Bilateral parameter (number of iterations)"},
    {"cleanup",VBooleanRepn,1,(VPointer) &cleanup,VOptionalOpt,NULL,"Whether to remove isloated voxels"},
    {"mask", VStringRepn, 1, (VPointer) &mask_filename, VRequiredOpt, NULL, "Mask"},
    {"j",VShortRepn,1,(VPointer) &nproc,VOptionalOpt,NULL,"number of processors to use, '0' to use all"},
  };

  FILE *fp=NULL;
  VString in_filename;
  VAttrList out_list=NULL,geolist=NULL;
  int i;
  char *prg_name=GetLipsiaName("vlisa_precoloring");
  fprintf (stderr, "%s\n", prg_name);

  /* parse command line */
  if (! VParseCommand (VNumber (options), options, & argc, argv)) {
    VReportUsage (argv[0], VNumber (options), options, NULL);
    exit (EXIT_FAILURE);
  }
  if (argc > 1) {
    VReportBadArgs (argc, argv);
    exit (EXIT_FAILURE);
  }

  /* omp-stuff */
#ifdef _OPENMP
  int num_procs=omp_get_num_procs();
  if (nproc > 0 && nproc < num_procs) num_procs = nproc;
  fprintf(stderr," using %d cores\n",(int)num_procs);
  omp_set_num_threads(num_procs);
#endif /* _OPENMP */

  /* read functional image data */
  int nlists = in_files.number;
  if (nlists < 1) VError(" no input");

  VAttrList *list = (VAttrList *) VCalloc(nlists,sizeof(VAttrList));
  for (i=0; i<nlists; i++) {
    in_filename = ((VString *) in_files.vector)[i];
    fprintf(stderr," %3d:  %s\n",i,in_filename);
    list[i] = VReadAttrList(in_filename,0L,TRUE,FALSE);
    if (geolist == NULL) {
      geolist = VGetGeoInfo(list[i]);
      double *DGeo = VGetGeoDim(geolist,NULL);
      if (fabs(DGeo[0]-4.0) > 0.01) VError(" Input files must be 4D (not 3D)");
    }
  }

  /* get number og design files */
  int dlists = des_files.number;
  if (dlists != nlists) {
    VError(" number of input functional files (%d) and design files (%d) do not match",nlists,dlists);
  }

  /* apply brain mask or threshold */
  VMultMinval(list,nlists,mask_filename,0.0);

  /* read data and voxel map */
  double tr=0;
  float *run_duration = (float *) VCalloc(nlists,sizeof(float));
  VGetTimeInfos(list,nlists,&tr,run_duration);
  gsl_matrix *Data = VReadImageData(list,nlists);
  VImage map  = VoxelMap(list[0]);
  int nslices = VPixel(map,0,3,0,VShort);
  int nrows   = VPixel(map,0,3,1,VShort);
  int ncols   = VPixel(map,0,3,2,VShort);
  int ntimesteps = Data->size2;

  /* additional regressors, no task labels, not included in permutations */
  gsl_matrix *ctmp1=NULL;
  gsl_matrix *ctmp2=NULL;
  gsl_matrix *covariates=NULL;
  int cdim = 1;
  int nuisance_dim=0;
  if (strlen(cova_filename) > 1) {
    ctmp1 = VReadCovariates(cova_filename,demean);
    if (ctmp1->size1 != Data->size2) VError(" num timesteps in covariate file not consistent with data");
    nuisance_dim = ctmp1->size2;
  }
  if (globalmean) {
    if (ctmp1 != NULL) cdim = ctmp1->size2+1;
    ctmp2 = gsl_matrix_calloc(Data->size2,cdim);
    GlobalMean(Data,ctmp2,(int)(cdim-1));
  }
  if (ctmp1 != NULL && ctmp2 == NULL) covariates = ctmp1;
  if (ctmp2 != NULL) covariates = ctmp2;

  /* design files with task labels */
  Trial **trial = (Trial **) VCalloc(dlists,sizeof(Trial *));
  int *numtrials = (int *) VCalloc(dlists,sizeof(int *));
  int nevents = 0;
  int sumtrials = 0;
  for (i=0; i<dlists; i++) {
    in_filename = ((VString *) des_files.vector)[i];
    fprintf(stderr," %3d:  %s\n",i,in_filename);
    int kk=0,jj=0;
    trial[i] = ReadDesign(in_filename,&kk,&jj);
    numtrials[i] = kk;
    if (jj > nevents) nevents = jj;
    sumtrials += numtrials[i];
  }
  fprintf(stderr," Number of trials: %d, number of event types: %d\n",sumtrials,nevents);
  Trial *alltrials = ConcatenateTrials(trial,numtrials,run_duration,nlists,sumtrials);
  CheckTrialLabels(alltrials,sumtrials);

  /* read contrast vector */
  gsl_vector *cont = gsl_vector_alloc(contrast.number + nuisance_dim);
  gsl_vector_set_zero(cont);
  for (i=0; i < contrast.number; i++) {
    double u = ((VFloat *)contrast.vector)[i];
    gsl_vector_set(cont,i,u);
  }

  /* alloc initial design matrix X */
  gsl_matrix *X = VCreateDesign(ntimesteps,nevents,(int)hemomodel,firstcol,covariates);
  fprintf(stderr," Design file dimensions: %lu x %lu\n",X->size1,X->size2);

  gsl_matrix *XInv = gsl_matrix_calloc(X->size2,X->size1);
  if (X->size2 != cont->size) {
    VError(" dimension of contrast vector (%ld) does not match design matrix (%ld)",cont->size-1,X->size2);
  }

  /* ini random permutations */
  gsl_rng_env_setup();
  const gsl_rng_type *T = gsl_rng_default;
  gsl_rng *rx = gsl_rng_alloc(T);
  gsl_rng_set(rx,(unsigned long int)seed);
  if (verbose) fprintf(stderr," seed: %ld\n",(long)seed);
  int **permtable = genperm(rx,numtrials,sumtrials,dlists,(int)numperm);

  /* estimate null variance to adjust radiometric parameter, use first 30 permutations */
  int nperm=0;
  float stddev = 1.0;
  double meanvar = 0.0;
  if (numperm > 0) {
    int tstperm = 30;
    if (tstperm > numperm) tstperm = numperm;
    VImage zmap = VCreateImage(nslices,nrows,ncols,VFloatRepn);
    double varsum=0,nx=0;
    for (nperm = 0; nperm < tstperm; nperm++) {
      Trial *permtrials = CopyTrials(alltrials,sumtrials);
      int j=0;
      for (j=0; j<sumtrials; j++) {
	int j0 = permtable[nperm][j];
	permtrials[j].id = alltrials[j0].id;
      }
      gsl_matrix *X = VCreateDesign(ntimesteps,nevents,(int)hemomodel,firstcol,covariates);
      gsl_matrix *XInv = gsl_matrix_calloc(X->size2,X->size1);
      VHemoModel(permtrials,sumtrials,nevents,ntimesteps,tr,(int)hemomodel,firstcol,X,covariates);
      VGLM(Data,X,XInv,cont,map,zmap);
      varsum += VImageVar(zmap);
      nx++;

      gsl_matrix_free(X);
      gsl_matrix_free(XInv);
      VFree(permtrials);
    }
    meanvar = varsum/nx;
    stddev = (float)(sqrt(meanvar));  /* update stddev */
    VDestroyImage(zmap);
  }

  /* no permutation */
  VImage zmap1 = VCreateImage(nslices,nrows,ncols,VFloatRepn);
  VCopyImageAttrs (map,zmap1);
  VImage dst1 = VCreateImageLike (zmap1);
  VHemoModel(alltrials,sumtrials,nevents,ntimesteps,tr,(int)hemomodel,firstcol,X,covariates);
  if (strlen(plot_filename) > 0) PlotDesign(X,tr,plot_filename);
  VGLM(Data,X,XInv,cont,map,zmap1);

  if (numperm == 0) {
    double z = VImageVar(zmap1);
    stddev = sqrt(z); /* update stddev */
  }
  float mode=0;
  if (numperm > 0) VZScale(zmap1,mode,stddev);
  VBilateralFilter(zmap1,dst1,(int)radius,(double)rvar,(double)svar,(int)numiter);

  /* ini histograms */
  double hmin=0,hmax=0;
  VGetHistRange(dst1,&hmin,&hmax);
  size_t nbins = 20000;
  gsl_histogram *hist0 = gsl_histogram_alloc (nbins);
  gsl_histogram_set_ranges_uniform (hist0,hmin,hmax);
  gsl_histogram *histz = gsl_histogram_alloc (nbins);
  gsl_histogram_set_ranges_uniform (histz,hmin,hmax);
  HistoUpdate(dst1,histz);

  /* random permutations */
#pragma omp parallel for shared(Data) schedule(dynamic)
  for (nperm = 0; nperm < numperm; nperm++) {
    if (nperm%5 == 0) fprintf(stderr," perm  %4d  of  %d\r",nperm,(int)numperm);

    /* randomly shuffle trial labels */
    Trial *permtrials = CopyTrials(alltrials,sumtrials);
    int j=0;
    for (j=0; j<sumtrials; j++) {
      int j0 = permtable[nperm][j];
      permtrials[j].id = alltrials[j0].id;
    }

    /* hemodynamic model */
    gsl_matrix *X = VCreateDesign(ntimesteps,nevents,(int)hemomodel,firstcol,covariates);
    gsl_matrix *XInv = gsl_matrix_calloc(X->size2,X->size1);
    VHemoModel(permtrials,sumtrials,nevents,ntimesteps,tr,(int)hemomodel,firstcol,X,covariates);

    /* GLM */
    VImage zmap = VCreateImageLike(zmap1);
    VGLM(Data,X,XInv,cont,map,zmap);
    VZScale(zmap,mode,stddev);
    gsl_matrix_free(X);
    gsl_matrix_free(XInv);
    VFree(permtrials);

    /* bilateral filter */
    VImage dst = VCreateImageLike (zmap);
    VBilateralFilter(zmap,dst,(int)radius,(double)rvar,(double)svar,(int)numiter);

#pragma omp critical
    {
      HistoUpdate(dst,hist0);
    }
    VDestroyImage(dst);
    VDestroyImage(zmap);
  }

  /* apply fdr */
  VImage fdrimage = VCopyImage (dst1,NULL,VAllBands);
  if (numperm > 0) {
    FDR(dst1,fdrimage,hist0,histz,(double)alpha);
    if (cleanup && alpha < 1.0) {
      VIsolatedVoxels(fdrimage,(float)(1.0-alpha));
    }
  }

  /*
  ** output
  */
  out_list = VCreateAttrList ();
  VHistory(VNumber(options),options,prg_name,&list[0],&out_list);

  /* update geoinfo, 4D to 3D */
  if (geolist != NULL) {
    double *D = VGetGeoDim(geolist,NULL);
    D[0] = 3;
    D[4] = 1;
    VSetGeoDim(geolist,D);
  }
  VSetGeoInfo(geolist,out_list);
  VAppendAttr (out_list,"image",NULL,VImageRepn,fdrimage);
  fp = VOpenOutputFile (out_filename, TRUE);
  if (! VWriteFile (fp, out_list)) exit (1);
  fclose(fp);

  fprintf (stderr, "\n%s: done.\n", argv[0]);
  exit(0);
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
/* PSICOV - Protein Sparse Inverse COVariance analysis program */

/* by David T. Jones August 2011 - Copyright (C) 2011 University College London */

/* Version 1.05 - Last Edit 13/2/12 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <math.h>
#include <unistd.h>

#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>

#define FALSE 0
#define TRUE 1

#define SQR(x) ((x)*(x))
#define MAX(x,y) ((x)>(y)?(x):(y))
#define MIN(x,y) ((x)<(y)?(x):(y))

#define MAXSEQLEN 5000
#define MINSEQS 50
#define MINEFSEQS 100

extern glasso_(int *, double *, double *, int *, int *, int *, int *, double *, int *, double *, double *, int *, double *, int *);

/* Dump a rude message to standard error and exit */
void
                fail(char *errstr)
{
    fprintf(stderr, "\n*** %s\n\n", errstr);
    exit(-1);
}

/* Convert AA letter to numeric code (0-21) */
int
                aanum(int ch)
{
    const static int aacvs[] =
    {
	999, 0, 3, 4, 3, 6, 13, 7, 8, 9, 21, 11, 10, 12, 2,
	21, 14, 5, 1, 15, 16, 21, 19, 17, 21, 18, 6
    };

    return (isalpha(ch) ? aacvs[ch & 31] : 20);
}

/* Allocate matrix */
void           *allocmat(int rows, int columns, int size)
{
    int             i;
    void          **p, *rp;

    rp = malloc(rows * sizeof(void *) + sizeof(int));

    if (rp == NULL)
	fail("allocmat: malloc [] failed!");

    *((int *)rp) = rows;

    p = rp + sizeof(int);

    for (i = 0; i < rows; i++)
	if ((p[i] = calloc(columns, size)) == NULL)
	    fail("allocmat: malloc [][] failed!");

    return p;
}

/* Free matrix */
void
                freemat(void *rp)
{
    int             rows;
    void **p = rp;

    rows = *((int *)(rp - sizeof(int)));

    while (rows--)
	free(p[rows]);

    free(rp - sizeof(int));
}

/* Allocate vector */
void           *allocvec(int columns, int size)
{
    void          *p;

    p = calloc(columns, size);

    if (p == NULL)
	fail("allocvec: calloc failed!");

    return p;
}

struct sc_entry
{
    float sc;
    int i, j;
} *sclist;

/* Sort descending */
int cmpfn(const void *a, const void *b)
{
    if (((struct sc_entry *)a)->sc == ((struct sc_entry *)b)->sc)
	return 0;

    if (((struct sc_entry *)a)->sc < ((struct sc_entry *)b)->sc)
	return 1;

    return -1;
}

int             main(int argc, char **argv)
{
    int             a, b, i, j, k, seqlen, nids, s, nseqs, ncon, opt, ndim, approxflg=0, initflg=0, debugflg=0, diagpenflg=1, apcflg=1, maxit=10000, npair, nnzero, niter, jerr, shrinkflg=1, rawscflg = 1, pseudoc = 1, minseqsep = 5;
    unsigned int *wtcount;
    double thresh=1e-4, del, sum, score, (**pab)[21][21], **pa, wtsum, pc, **pcmat, *pcsum, pcmean, rhodefault = -1.0, lambda, smean, fnzero, lastfnzero, trialrho, rfact, r2, targfnzero = 0.0, scsum, scsumsq, mean, sd, zscore, ppv;
    float *weight, idthresh = -1.0, maxgapf = 0.9;
    char            buf[4096], seq[MAXSEQLEN], *blockfn = NULL, **aln;
    FILE *ifp;

    while ((opt = getopt(argc, argv, "alnpr:b:i:t:c:g:d:j:")) >= 0)
	switch (opt)
	{
	case 'a':
	    approxflg = 1;
	    break;
	case 'n':
	    shrinkflg = 0;
	    break;
	case 'p':
	    rawscflg = 0;
	    break;
	case 'l':
	    apcflg = 0;
	    break;
	case 'r':
	    rhodefault = atof(optarg);
	    break;
	case 'd':
	    targfnzero = atof(optarg);
	    break;
	case 't':
	    thresh = atof(optarg);
	    break;
	case 'i':
	    idthresh = 1.0 - atof(optarg)/100.0;
	    break;
	case 'c':
	    pseudoc = atoi(optarg);
	    break;
	case 'j':
	    minseqsep = atoi(optarg);
	    break;
	case 'b':
	    blockfn = strdup(optarg);
	    break;
	case 'g':
	    maxgapf = atof(optarg);
	    break;
	case '?':
	    exit(-1);
	}

    if (optind >= argc)
	fail("Usage: psicov [options] alnfile\n\nOptions:\n-a\t: use approximate Lasso algorithm\n-n\t: don't pre-shrink the sample covariance matrix\n-p\t: output PPV estimates rather than raw scores\n-l\t: don't apply APC to Lasso output\n-r nnn\t: set initial rho parameter\n-d nnn\t: set target precision matrix sparsity (default 0 = not specified)\n-t nnn\t: set Lasso convergence threshold (default 1e-4)\n-i nnn\t: select BLOSUM weighting with given identity threshold (default selects threshold automatically)\n-c nnn\t: set pseudocount value (default 1)\n-j nnn\t: set minimum sequence separation (default 5)\n-g nnn\t: maximum fraction of gaps (default 0.9)\n-b file\t: read rho parameter file\n");

    ifp = fopen(argv[optind], "r");
    if (!ifp)
	fail("Unable to open alignment file!");

    for (nseqs=0;; nseqs++)
	if (!fgets(seq, MAXSEQLEN, ifp))
	    break;

    aln = allocvec(nseqs, sizeof(char *));

    weight = allocvec(nseqs, sizeof(float));

    wtcount = allocvec(nseqs, sizeof(unsigned int));

    rewind(ifp);

    if (!fgets(seq, MAXSEQLEN, ifp))
	fail("Bad alignment file!");

    seqlen = strlen(seq)-1;

    if (nseqs < MINSEQS)
	fail("Alignment too small - not enough homologous sequences to proceed (or change MINSEQS at your own risk!)");

    if (!(aln[0] = malloc(seqlen)))
	fail("Out of memory!");

    for (j=0; j<seqlen; j++)
	aln[0][j] = aanum(seq[j]);

    for (i=1; i<nseqs; i++)
    {
	if (!fgets(seq, MAXSEQLEN, ifp))
	    break;

	if (seqlen != strlen(seq)-1)
	    fail("Length mismatch in alignment file!");

	if (!(aln[i] = malloc(seqlen)))
	    fail("Out of memory!");

	for (j=0; j<seqlen; j++)
	    aln[i][j] = aanum(seq[j]);
    }

    /* Calculate sequence weights */

    if (idthresh < 0.0)
    {
	double meanfracid = 0.0;

	for (i=0; i<nseqs; i++)
	    for (j=i+1; j<nseqs; j++)
	    {
		int nids;
		float fracid;

		for (nids=k=0; k<seqlen; k++)
		    if (aln[i][k] == aln[j][k])
			nids++;

		fracid = (float)nids / seqlen;

		meanfracid += fracid;
	    }

	meanfracid /= 0.5 * nseqs * (nseqs - 1.0);

	idthresh = 0.38 * 0.32 / meanfracid;
    }

    for (i=0; i<nseqs; i++)
	for (j=i+1; j<nseqs; j++)
	{
	    int nthresh = seqlen * idthresh;

	    for (k=0; nthresh > 0 && k<seqlen; k++)
		if (aln[i][k] != aln[j][k])
		    nthresh--;

	    if (nthresh > 0)
	    {
		wtcount[i]++;
		wtcount[j]++;
	    }
	}

    for (wtsum=i=0; i<nseqs; i++)
	wtsum += (weight[i] = 1.0 / (1 + wtcount[i]));

    if (wtsum < MINEFSEQS)
	puts("\n*** WARNING - not enough sequence variation - or change MINEFSEQS at your own risk! ***\n");

    pa = allocmat(seqlen, 21, sizeof(double));
    pab = allocmat(seqlen, seqlen, 21*21*sizeof(double));

    /* Calculate singlet frequencies with pseudocount */
    for (i=0; i<seqlen; i++)
    {
	for (a=0; a<21; a++)
	    pa[i][a] = pseudoc;

	for (k=0; k<nseqs; k++)
	{
	    a = aln[k][i];
	    if (a < 21)
		pa[i][a] += weight[k];
	}

	for (a=0; a<21; a++)
	    pa[i][a] /= pseudoc * 21.0 + wtsum;
    }

    /* Calculate pair frequencies with pseudocount */
    for (i=0; i<seqlen; i++)
    {
	for (j=i+1; j<seqlen; j++)
	{
	    for (a=0; a<21; a++)
		for (b=0; b<21; b++)
		    pab[i][j][a][b] = pseudoc / 21.0;

	    for (k=0; k<nseqs; k++)
	    {
		a = aln[k][i];
		b = aln[k][j];
		if (a < 21 && b < 21)
		    pab[i][j][a][b] += weight[k];
	    }

	    for (a=0; a<21; a++)
		for (b=0; b<21; b++)
		{
		    pab[i][j][a][b] /= pseudoc * 21.0 + wtsum;
		    pab[j][i][b][a] = pab[i][j][a][b];

//		    printf("%d/%d %d/%d %f %f %f %f\n", i+1, a, j+1, b, pab[i][j][a][b], pa[i][a] , pa[j][b], pab[i][j][a][b] - pa[i][a] * pa[j][b]);
		}
	}
    }

    for (i=0; i<seqlen; i++)
	for (a=0; a<21; a++)
	    for (b=0; b<21; b++)
		pab[i][i][a][b] = (a == b) ? pa[i][a] : 0.0;

    gsl_matrix *cmat, *rho, *ww, *wwi, *tempmat;

    ndim = seqlen * 21;

    cmat = gsl_matrix_calloc(ndim, ndim);

    /* Form the covariance matrix */
    for (i=0; i<seqlen; i++)
	for (j=0; j<seqlen; j++)
	    for (a=0; a<21; a++)
		for (b=0; b<21; b++)
		    if (i != j)
			gsl_matrix_set(cmat, i*21+a, j*21+b, pab[i][j][a][b] - pa[i][a] * pa[j][b]);
		    else if (a == b)
			gsl_matrix_set(cmat, i*21+a, j*21+b, pab[i][j][a][b] - pa[i][a] * pa[j][b]);

    freemat(pab);

    /* Shrink sample covariance matrix towards shrinkage target F = Diag(1,1,1,...,1) * smean */

    if (shrinkflg)
    {
	for (smean=i=0; i<ndim; i++)
	    smean += gsl_matrix_get(cmat, i, i);

	smean /= (float)ndim;
	lambda = 0.1;

//	smean = 1;

	tempmat = gsl_matrix_calloc(ndim, ndim);

	gsl_set_error_handler_off();

	for (;;)
	{
	    gsl_matrix_memcpy(tempmat, cmat);

	    /* Test if positive definite using Cholesky decomposition */
	    if (!gsl_linalg_cholesky_decomp(tempmat))
		break;

	    for (i=0; i<seqlen; i++)
		for (j=0; j<seqlen; j++)
		    for (a=0; a<21; a++)
			for (b=0; b<21; b++)
			    if (i != j)
				gsl_matrix_set(cmat, i*21+a, j*21+b, (1.0 - lambda) * gsl_matrix_get(cmat, i*21+a, j*21+b));
			    else if (a == b)
				gsl_matrix_set(cmat, i*21+a, j*21+b, smean * lambda + (1.0 - lambda) * gsl_matrix_get(cmat, i*21+a, j*21+b));
	}

	gsl_matrix_free(tempmat);
    }

    rho = gsl_matrix_alloc(ndim, ndim);
    ww = gsl_matrix_alloc(ndim, ndim);
    wwi = gsl_matrix_alloc(ndim, ndim);

    lastfnzero=0.0;

    /* Guess at a reasonable starting rho value if undefined */
    if (rhodefault < 0.0)
	trialrho = MAX(0.001, 1.0 / wtsum);
    else
	trialrho = rhodefault;

    rfact = 0.0;

    for (;;)
    {
	if (trialrho <= 0.0 || trialrho >= 1.0)
	    fail("Sorry - failed to find suitable value for rho (0 < rho < 1)!");

	gsl_matrix_set_all(rho, trialrho);

	for (i=0; i<seqlen; i++)
	    for (j=0; j<seqlen; j++)
		for (a=0; a<21; a++)
		    for (b=0; b<21; b++)
			if ((a != b && i == j) || pa[i][20] > maxgapf || pa[j][20] > maxgapf)
			    gsl_matrix_set(rho, i*21+a, j*21+b, 1e9);

	/* Mask out regions if block-out list provided */
	if (blockfn != NULL)
	{
	    ifp = fopen(blockfn, "r");

	    for (;;)
	    {
		if (fscanf(ifp, "%d %d %lf", &i, &j, &score) != 3)
		    break;

		for (a=0; a<21; a++)
		    for (b=0; b<21; b++)
		    {
			gsl_matrix_set(rho, (i-1)*21+a, (j-1)*21+b, score);
			gsl_matrix_set(rho, (j-1)*21+b, (i-1)*21+a, score);
		    }
	    }

	    fclose(ifp);
	}

	/* All matrices are symmetric so no need to transpose before/after calling Fortran code */

	glasso_(&ndim, cmat->data, rho->data, &approxflg, &initflg, &debugflg, &diagpenflg, &thresh, &maxit, ww->data, wwi->data, &niter, &del, &jerr);

	if (targfnzero <= 0.0)
	    break;

	for (npair=nnzero=i=0; i<ndim; i++)
	    for (j=i+1; j<ndim; j++,npair++)
		if (gsl_matrix_get(wwi, i, j) != 0.0)
		    nnzero++;

	fnzero = (double) nnzero / npair;

//      printf("rho=%f fnzero = %f\n", trialrho, fnzero);

	/* Stop iterating if we have achieved the target sparsity level */
	if (fabs(fnzero - targfnzero)/targfnzero < 0.01)
	    break;

	if (fnzero == 0.0)
	{
	    /* As we have guessed far too high, halve rho and try again */
	    trialrho *= 0.5;
	    continue;
	}

	if (lastfnzero > 0.0 && fnzero != lastfnzero)
	{
//	    printf("fnzero=%f lastfnzero=%f trialrho=%f oldtrialrho=%f\n", fnzero, lastfnzero, trialrho, trialrho/rfact);

	    rfact = pow(rfact, log(targfnzero / fnzero) / log(fnzero / lastfnzero));

//	    printf("New rfact = %f\n", rfact);
	}

	lastfnzero = fnzero;

	/* Make a small trial step in the appropriate direction */

	if (rfact == 0.0)
	    rfact = (fnzero < targfnzero) ? 0.9 : 1.1;

	trialrho *= rfact;
    }

    gsl_matrix_free(rho);
    gsl_matrix_free(ww);

    /* Calculate background corrected scores using average product correction */

    pcmat = allocmat(seqlen, seqlen, sizeof(double));
    pcsum = allocvec(seqlen, sizeof(double));

    pcmean = 0.0;

    for (i=0; i<seqlen; i++)
	for (j=i+1; j<seqlen; j++)
	{
	    for (pc=a=0; a<20; a++)
		for (b=0; b<20; b++)
		    pc += fabs(gsl_matrix_get(wwi, i*21+a, j*21+b));

	    pcmat[i][j] = pcmat[j][i] = pc;
	    pcsum[i] += pc;
	    pcsum[j] += pc;

	    pcmean += pc;
	}

    pcmean /= seqlen * (seqlen - 1) * 0.5;

    /* Build final list of predicted contacts */

    sclist = allocvec(seqlen * (seqlen - 1) / 2, sizeof(struct sc_entry));

    for (scsum=scsumsq=ncon=i=0; i<seqlen; i++)
	for (j=i+minseqsep; j<seqlen; j++)
	    if (pcmat[i][j] > 0.0)
	    {
		/* Calculate APC score */
		if (apcflg)
		    sclist[ncon].sc = pcmat[i][j] - pcsum[i] * pcsum[j] / SQR(seqlen - 1.0) / pcmean;
		else
		    sclist[ncon].sc = pcmat[i][j];
		scsum += sclist[ncon].sc;
		scsumsq += SQR(sclist[ncon].sc);
		sclist[ncon].i = i;
		sclist[ncon++].j = j;
	    }

    qsort(sclist, ncon, sizeof(struct sc_entry), cmpfn);

    mean = scsum / ncon;
    sd = sqrt(scsumsq / ncon - SQR(mean));

    /* Print output in CASP RR format with optional PPV estimated from final Z-score */
    if (rawscflg)
	for (i=0; i<ncon; i++)
	    printf("%d %d 0 8 %f\n", sclist[i].i+1, sclist[i].j+1, sclist[i].sc);
    else
	for (i=0; i<ncon; i++)
	{
	    zscore = (sclist[i].sc - mean) / sd;
	    ppv = 0.904 / (1.0 + 16.61 * exp(-0.8105 * zscore));
	    printf("%d %d 0 8 %f\n", sclist[i].i+1, sclist[i].j+1, ppv);
	}

    return 0;
}


text:
#include <lapacke.h>

#include "tasks.h"

void chol_task_seq(void *ptr)
{
    struct chol_task_arg *arg = (struct chol_task_arg*) ptr;

    int n     = arg->n;
    double *A = arg->A;
    int ldA   = arg->ldA;

    LAPACKE_dpotrf(LAPACK_COL_MAJOR, 'L', n, A, ldA);
}

max_stars_repo_name: NLAFET/pcp-runtime

max_issues_repo_name: NLAFET/pcp-runtime

max_forks_repo_name: NLAFET/pcp-runtime


text:
/* Copyright (c) 2011-2012, Jérémy Fix. All rights reserved. */

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

/* * Redistributions of source code must retain the above copyright notice, */
/* this list of conditions and the following disclaimer. */
/* * Redistributions in binary form must reproduce the above copyright notice, */
/* this list of conditions and the following disclaimer in the documentation */
/* and/or other materials provided with the distribution. */
/* * None of the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. */

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

#ifndef EKF_TYPES_H
#define EKF_TYPES_H

#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include "ukf_math.h"

/**
 * @short Extended Kalman Filter in the case of additive noise, the notations follow Van Der Merwe, phD, p. 36
 */
namespace ekf
{

  class EvolutionNoise;

  /**
   * @short Structure holding the parameters of the statistical linearization
   *
   */
  typedef struct
  {

    /**
     * @short Evolution noise type
     */
    EvolutionNoise *evolution_noise;

    /**
     * @short Covariance of the observation noise
     */
    double observation_noise;

    /**
     * @short Prior estimate of the covariance matrix
     */
    double prior_pk;

    /**
     * @short Number of parameters to estimate
     */
    int n;

    /**
     * @short Dimension of the output
     */
    int no;

    /**
     * @short Is the observation gradient diagonal ? In that case, simplifications can be introduced
     */
    bool observation_gradient_is_diagonal;

  } ekf_param;

  typedef struct
  {
    /**
     * @short Current estimate of the state
     */
    gsl_vector * xk;

    /**
     * @short Predicted state
     */
    gsl_vector * xkm;

    /**
     * @short Variance-Covariance of the state
     */
    gsl_matrix * Pxk;

    /**
     * @short Jacobian of the evolution function
     */
    gsl_matrix *Fxk;

    /**
     * @short Variance covariance of the evolution noise
     */
    gsl_matrix * Rv;

    /**
     * @short Current observation
     */
    gsl_vector * yk;

    /**
     * @short Current innovations
     */
    gsl_vector * ino_yk;

    /**
     * @short Jacobian of the observation function
     */
    gsl_matrix *Hyk;

    /**
     * @short Variance covariance of the observation noise
     */
    gsl_matrix * Rn;

    /**
     * @short Kalman gain
     */
    gsl_matrix *Kk;

    /**
     * @short Temporary matrices
     */
    gsl_matrix * temp_n_n;
    gsl_matrix * temp_n_1;
    gsl_matrix * temp_no_no;
    gsl_matrix * temp_n_no;
    gsl_matrix * temp_2_n_n;
    gsl_vector * temp_no;

    /**
     * @short Optional parameters (this must be allocated and initialized from the user side!
     */
    gsl_vector * params;

  } ekf_state;

  /**
   * @short Mother class from which the evolution noises inherit
   */
  class EvolutionNoise
  {
  protected:
    double _initial_value;
  public:
  EvolutionNoise(double initial_value) : _initial_value(initial_value) {};
    void init(ekf_param &p, ekf_state &s)
    {
      gsl_matrix_set_identity(s.Rv);
      gsl_matrix_scale(s.Rv, _initial_value);
    }

    virtual void updateEvolutionNoise(ekf_param &p, ekf_state &s) = 0;
  };

  /**
   * @short Annealing type evolution noise
   */
  class EvolutionAnneal : public EvolutionNoise
  {
    double _decay, _lower_bound;
  public:
  EvolutionAnneal(double initial_value, double decay, double lower_bound) : EvolutionNoise(initial_value),
      _decay(decay),
      _lower_bound(lower_bound)
      { };

    void updateEvolutionNoise(ekf_param &p, ekf_state &s)
    {
      for(int i = 0 ; i < p.n ; ++i)
	{
	  for(int j = 0 ; j < p.n ; ++j)
	    {
	      if(i == j)
		gsl_matrix_set(s.Rv, i, i, ukf::math::max(_decay * gsl_matrix_get(s.Rv,i,i),_lower_bound));
	      else
		gsl_matrix_set(s.Rv, i, j, 0.0);
	    }
	}
    }
  };

  /**
   * @short Forgetting type evolution noise
   */
  class EvolutionRLS : public EvolutionNoise
  {
    double _decay;
  public:
  EvolutionRLS(double initial_value, double decay) : EvolutionNoise(initial_value), _decay(decay)
    {
      if(ukf::math::cmp_equal(_decay, 0.0))
	printf("Forgetting factor should not be null !!\n");
    };

    void updateEvolutionNoise(ekf_param &p, ekf_state &s)
    {
      gsl_matrix_memcpy(s.Rv, s.Pxk);
      gsl_matrix_scale(s.Rv, 1.0 / _decay - 1.0);
    }
  };

  /**
   * @short Robbins-Monro evolution noise
   */
  class EvolutionRobbinsMonro : public EvolutionNoise
  {
    double _alpha;
  public:
  EvolutionRobbinsMonro(double initial_value, double alpha) : EvolutionNoise(initial_value),
      _alpha(alpha)
      { };

    void updateEvolutionNoise(ekf_param &p, ekf_state &s)
    {
      // Compute Kk * ino_yk
      gsl_matrix_view mat_view = gsl_matrix_view_array(s.ino_yk->data, p.no, 1);
      gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, s.Kk, &mat_view.matrix, 0.0, s.temp_n_1);
      // Compute : Rv = (1 - alpha) Rv + alpha . (Kk * ino_yk) * (Kk * ino_yk)^T
      gsl_blas_dgemm(CblasNoTrans, CblasTrans, _alpha, s.temp_n_1, s.temp_n_1, (1.0 - _alpha),s.Rv);
    }
  };

}

#endif // EKF_TYPES_H


text:
#ifndef CWANNIER_PARTIALNUMVALUES_H
#define CWANNIER_PARTIALNUMVALUES_H

#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_eigen.h>
#include "HTightBinding.h"
#include "PartialDosValues.h"
#include "ctetra/sum.h"

double** PartialNumValues(HTightBinding *Hrs, gsl_matrix *R, int na, int nb, int nc, double num_total_electrons, double **Es, double num_E, double *E_Fermi, double **num_states_Fermi);

#endif //CWANNIER_PARTIALNUMVALUES_H

max_stars_repo_name: tflovorn/cwannier

max_issues_repo_name: tflovorn/cwannier

max_forks_repo_name: tflovorn/cwannier


text:
/*
F(4,3)
cc -static -o test_cblas_open main43.c -I /opt/OpenBLAS/include/ -L/opt/OpenBLAS/lib -lopenblas -lpthread -lgfortran
	->For each cube of kernel
		-> for each tile (a tile is a cube d x d x channel)
			-> for each channel
				-> Apply the winograd algorithm
			-> sum the result of the previous channel to the result of the next channell
		-> write the tile in a single output layer
	->Add the output layer in che cube layer and change the kernel cube
*/

#include <cblas.h>
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <time.h>
#include <string.h>

#define MAX_ELEMENTS 100
#define FILENAMELEN 10
#define DDIMENSION(m,r) (m+r-1)
#define OUTPUTDIMENSION(m,w,d) (m * (w / d))

typedef float * Matrix;
typedef Matrix * Cube;
typedef Cube * Hypercube;

void productABAt(float result[], float A[], float B[], int rowsA, int colsB, int colsA);
void twoDcorrelation(float C[], float g[], float D[], int m, int r, int d);
Matrix elaborateTile(Cube tile, Cube g, int m, int r, int d, int channel);
Cube fillTheTile(Cube image, int starth, int startw, int w, int d, int channel);
void saveOutputTile(Matrix output, Matrix tileOutput, int starth, int startw, int w, int d, int m);
Matrix elaborateKernel(Cube image, Cube g, int m, int r, int channel, int w, int h, int d);

/*START FUNCTION TO CREATE VARIABLES*/
	Cube generateCube(int w, int h, int ch);
	Matrix generateMatrix(int w, int h);
	Hypercube readKernels(int w, int h, int ch, int k, char fileName[]);
	Cube readInput(int w, int h, int ch, char fileName[]);

	void printCube(Cube c, int w, int channel);
/*END FUNCTION TO CREARE VARIABLES*/

/*START FUNCTIONS TO READ THE FILES*/
	void readFile(float result[], char fileName[]);
	void generateNameAT(char* at, int m, int r);
	void generateNameG(char* g, int m, int r);
	void generateNameBT(char* g, int m, int r);
/* END FUNCRIONS TO READ THE FILES*/

static Matrix A;	//Parameters that are used to make the lowest operation.
static Matrix B;	//
static Matrix G;	//

void main(int argc, char *argv[])
{

	int m = atoi(argv[1]);			//dimension of the output tile
	int r = atoi(argv[2]);			//dimension of the kernel rxr
	int channel = atoi(argv[3]);	//dimension of the channel
	int k = atoi(argv[4]);			//number of kernel
	int w = atoi(argv[5]);			//width of the image -> multiple of d
	int h = atoi(argv[6]);			//height of the image -> multiple of d

	char * inputFilename = argv[7];
	char * kernelFilename = argv[8];

	int d = DDIMENSION(m,r);

	char fnameAT[FILENAMELEN];
	char fnameG[FILENAMELEN];
	char fnameBT[FILENAMELEN];

	generateNameAT(fnameAT, m,r);
	generateNameG(fnameG, m,r);
	generateNameBT(fnameBT, m,r);

	A = generateMatrix(m, d);
	B = generateMatrix(d, d);
	G = generateMatrix(r, d);

	readFile(A, fnameAT);
	readFile(B, fnameBT);
	readFile(G, fnameG);

	Cube image = readInput(w, h, channel, inputFilename);
	Hypercube kernel = readKernels(r, r, channel, k, kernelFilename);

	printf("IMAGE\n");
	printCube(image, w, channel);

	printf("\nKERNELS\n");
	for(int i = 0; i < k; i++)
	{
		printf("Kernel number %d\n", i);
		printCube(kernel[i], r, channel);
	}
	printf("\n------------------\n");
	Cube output = (Cube)malloc(k * sizeof(Matrix));	//k is the number of filters

	//printf("START THE COMPUTATION - SLEEP FOR 3 SECONDS\n");
 	//sleep(3);
	/*
	 *	START COMPUTATION.
	 */
	clock_t start = clock();
	for(int i = 0; i < k; i ++)
		output[i] = elaborateKernel(image, kernel[i], m, r, channel, w, h, d);

	clock_t end = clock();
	/*
	 *	END COMPUTATION.
	 */

	printf("OUTPUT:\n");
	printCube(output, OUTPUTDIMENSION(m,w,d), k);

	float seconds = (float)(end - start) / CLOCKS_PER_SEC;

	printf("\nTIME IN SECONDS TO DO F(%d,%d) = %f\n",m,r,seconds);

}

void printCube(Cube c, int w, int channel)
{
	for(int i = 0; i < channel; i++)
	{
		printf("channel %d:\n",i);
		for(int j = 0; j < w * w; j ++)
		{
			printf("%.2f ", c[i][j]);
			if((j+1) % w == 0) printf("\n");
		}
	}
}

Matrix elaborateKernel(Cube image, Cube g, int m, int r, int channel, int w, int h, int d)
{
	Matrix output = generateMatrix(OUTPUTDIMENSION(m,w,d), OUTPUTDIMENSION(m,w,d));

	for(int i = 0; i < h; i = i + d)
	{
		for (int j = 0; j < w; j = j + d)
		{
			Cube tile = fillTheTile(image, i, j, w, d, channel);

			Matrix tileOutput = elaborateTile(tile, g, m, r, d, channel);
			saveOutputTile(output, tileOutput, i, j, w, d, m);

			free(tileOutput);
			free(tile);
		}
	}
	return output;
}

/*
	Write the output tile on the output image.
*/

void saveOutputTile(Matrix output, Matrix tileOutput, int starth, int startw, int w, int d, int m)
{
		int outputWidth = OUTPUTDIMENSION(m,w,d);
		int offsetw = (startw / d) * m;
		int offseth = (starth / d) * m;
		for(int y = 0; y < m; y ++)
			for(int x = 0; x < m; x++)
				output[(x + offsetw) + (outputWidth * (y + offseth))] = tileOutput[x + (m * y)];

}

/*
	This function takes an image, it divides the image based on starth and startw and it returns a tile.
		- image  -> the all image
		- starth -> height coordinate of the image where to start to make the tile
		- startw -> width coordinate of the image where to start to make the tile
		- w -> total width of the image
		- d -> dimension of the tile -> d x d
		- channel -> number of channels of the image
	RESULT
		- A cube d x d x channell that represent a single tile.
*/
Cube fillTheTile(Cube image, int starth, int startw, int w, int d, int channel)
{
	Cube tile = generateCube(d, d, channel);
	for(int c = 0; c < channel; c++)
		for(int y = 0; y < d; y ++)
			for(int x = 0; x < d; x++)
				tile[c][x + (d * y)] = image[c][(startw + x) + (w * (y + starth))];
	return tile;
}

/*
	This function is used to make the combination with a tile and a kernel.
		- tile has the dimension d x d and has 'channel' channels
		- g is the kernel, it has the dimension of r x r and has 'channel' channels
		- m -> dimension of the output
		- r -> dimension of the kernel
		- d -> dimension of the tile
		- channell -> number of channels
	RESULT
		- it returns a matrix calculated by an element wise sum of the result of the function twoDcorrelation
			calculated for each channel.
*/
Matrix elaborateTile(Cube tile, Cube g, int m, int r, int d, int channel)
{
	Matrix temp = generateMatrix(d,d);
	Matrix tileOutputBig = generateMatrix(d,d);
	Matrix tileOutput = generateMatrix(m,m);

	for(int i = 0; i < m * m; i++)
			tileOutput[i] = 0;

	for(int i = 0; i < channel; i++)
	{
	 	twoDcorrelation(temp, g[i], tile[i], m, r, d);
		for(int j = 0 ; j < d * d; j++)
			tileOutputBig[j] = tileOutputBig[j] + temp[j];

	}
	productABAt(tileOutput, A, tileOutputBig, m, d, d);
	free(temp);
	free(tileOutputBig);

	return tileOutput;
}

Cube generateCube(int w, int h, int ch)
{
	Cube cube = (Cube)malloc(ch * sizeof(Matrix));
	for(int i = 0; i < ch; i++)
		cube[i] = generateMatrix(w,h);

	return cube;
}

Matrix generateMatrix(int w, int h)
{
	return (Matrix)malloc(w * h * sizeof(float));
}

/*
	Calculate a single channel for a single kernel for a single tile
	Input Parameter:
		- result -> matrix of the result
		- g -> the considered kernel
		- D -> the considered tile
		- m -> Dimension of the output
		- r -> dimension of the kernel
		- d -> dimension of the tile (m + r - 1)
*/
void twoDcorrelation(float C[], float g[], float D[], int m, int r, int d)
{

	float GgG[MAX_ELEMENTS];
	float BDB[MAX_ELEMENTS];
	int i=0;

	productABAt(GgG, G, g, d, r, r);
	productABAt(BDB, B, D, d, d, d);

 	for(i=0;i<d*d;i++)
		C[i] = GgG[i] * BDB[i];
  	//productABAt(C, A, GgG, m, d, d);
}

/*
 *   Make the operation R = ABA**T
 *   Parameters:
 *		- Matrix Result
 *		- Matrix A and B
 *		- Number of rows of A -> rowsA
 *		- Number of cols of B -> colsB
 *		- Number of cols of A that is equal to the numer of rows of B -> colsA
 *  Result
 *		- Matrix result -> that is an rowsA x rowsA matrix
 */
void productABAt(float result[], float A[], float B[], int rowsA, int colsB, int colsA)
{
	float C[MAX_ELEMENTS];

	cblas_sgemm(CblasRowMajor,	//Modo in cui è salvata la matrice, legge i numeri riga per riga
		CblasNoTrans,	//Non fare trasposta
		CblasNoTrans, 	//Non fare trasposta
		rowsA, 		//numero righe matrice A
		colsB, 		//numero colonne matrice B
		colsA, 		//numero colonne A & numero righe B
		1, 		//moltiplicatore prima matrice
		A, 		//Matrice A
		colsA, 		//numero colonne matrice A
		B, 		//Matrice B
		colsB, 		//numero colonne matrice B
		0, 		//Moltiplicatore matrice C
		C,		//matrice C
		colsB);		//numero righe matrice C

	//C has rowsA rows and colsB cols
	//the result has rowsA rows and rowsA cols

	cblas_sgemm(CblasRowMajor,
		CblasNoTrans,
		CblasTrans,
		rowsA,	//numero righe matrice C
		rowsA,
		colsB,
		1,
		C,
		colsB,
		A,
		colsA,
		0,
		result,
		rowsA);
}

Cube readInput(int w, int h, int ch, char fileName[])
{

    FILE *f;
    f = fopen(fileName, "r");
    if(f==NULL)
	{
    	printf("I can not read the file %s\n", fileName);
		exit(0);
	}

    int fi = 0;
    int linesize = w*h;
    Cube matrix = (Cube)malloc(ch * sizeof(Matrix));

    for(int i = 0; i < ch; i++)
	{
        matrix[i] = (Matrix)malloc(linesize * sizeof(float));
        while(!feof(f)&&fi<linesize)
		{
            fscanf(f, "%f", &matrix[i][fi]);
            fi++;
        }
        fi=0;
    }
    fclose(f);
    return matrix;

}

Hypercube readKernels(int w, int h, int ch, int k, char fileName[])
{

    FILE *f;
    f = fopen(fileName, "r");
    if(f==NULL)
    {
		printf("I can not read the file %s\n", fileName);
		exit(0);
	}
    int fi = 0;
    int linesize = w*h;
    int kernelsize = linesize*ch;
    Hypercube matrix = (Hypercube)malloc(k * sizeof(Cube));

    for(int i = 0; i < k; i++)
	{
        matrix[i] = (Cube)malloc(kernelsize * sizeof(Matrix));
        for(int j=0; j< ch; j++)
		{
            matrix[i][j] = (Matrix)malloc(linesize * sizeof(float));
            while(!feof(f) && fi<linesize)
			{
                   fscanf(f, "%f", &matrix[i][j][fi]);
                    fi++;
            }
            fi=0;
        }
    }
    fclose(f);
    return matrix;
}

void readFile(float result[], char fileName[])
{
	FILE *f;
	f = fopen(fileName, "r");
	if(f == NULL)
	{
		printf("I can not read the file %s\n", fileName);
		exit(0);
	}
	int i = 0;
	while(!feof(f))
	{
		fscanf(f, "%f", &result[i]);
		i++;
	}
	fclose(f);
}

/* FILENAMES FUNCTIONS */
void generateNameAT(char* at, int m, int r){
    at[0] = 'A';
    at[1] = 'T';
    at[2] = m + '0';
    at[3] = 'x';
    at[4] = r + '0';
    at[5] = '.';
    at[6] = 't';
    at[7] = 'x';
    at[8] = 't';
    at[9] = '\0';
}

void generateNameG(char* g, int m, int r){
    g[0] = 'G';
    g[1] = m + '0';
    g[2] = 'x';
    g[3] = r + '0';
    g[4] = '.';
    g[5] = 't';
    g[6] = 'x';
    g[7] = 't';
    g[8] = '\0';
}

void generateNameBT(char* g, int m, int r){
    g[0] = 'B';
    g[1] = 'T';
    g[2] = m + '0';
    g[3] = 'x';
    g[4] = r + '0';
    g[5] = '.';
    g[6] = 't';
    g[7] = 'x';
    g[8] = 't';
    g[9] = '\0';
}

max_stars_repo_name: whcjimmy/winograd-convolutional-nn

max_issues_repo_name: whcjimmy/winograd-convolutional-nn

max_forks_repo_name: whcjimmy/winograd-convolutional-nn


text:
/*
  NAME:
     proj_gauss_mixtures_IDL
  PURPOSE:
     run the projected gaussian mixtures algorithm from IDL (or python)
  CALLING SEQUENCE:
     see IDL wrapper
  INPUT:
     from IDL wrapper
  OUTPUT:
     updated model gaussians and average loglikelihood, see IDL WRAPPER
  REVISION HISTORY:
     2008-09-21 - Written Bovy
     2010-03-01 Added noproj option - Bovy
     2010-04-01 Added noweight option and logweights - Bovy
*/
#include <stdio.h>
#include <stdbool.h>
#include <time.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <proj_gauss_mixtures.h>

int proj_gauss_mixtures_IDL(double * ydata, double * ycovar,
			    double * projection, double * logweights,
			    int N, int dy,
			    double * amp, double * xmean,
			    double * xcovar, int d, int K,
			    char * fixamp, char * fixmean,
			    char * fixcovar,
			    double * avgloglikedata, double tol,
			    int maxiter, char likeonly, double w,
			    char * logfilename, int slen, int splitnmerge,
			    char * convlogfilename, int convloglen,
			    char noprojection,char diagerrors,
			    char noweights){
  //Set up logfiles
  bool keeplog = true;
  char logname[slen+1];
  char convlogname[convloglen+1];
  int ss;
  if (*logfilename == 0 || likeonly != 0 || slen == 0)
    keeplog = false;
  else {
    for (ss = 0; ss != slen; ++ss)
      logname[ss] = (char) *(logfilename++);
    for (ss = 0; ss != convloglen; ++ss)
      convlogname[ss] = (char) *(convlogfilename++);
    logfilename -= slen;
    convlogfilename -= convloglen;
    logname[slen] = '\0';
    convlogname[convloglen] = '\0';
  }

  if (keeplog) {
    logfile = fopen(logname,"a");
    if (logfile == NULL) return -1;
    convlogfile = fopen(convlogname,"w");
    if (convlogfile == NULL) return -1;
  }

  if (keeplog){
    time_t now;
    time(&now);
    fprintf(logfile,"#----------------------------------\n");
    fprintf(logfile,"#\n#%s\n",asctime(localtime(&now)));
    fprintf(logfile,"#----------------------------------\n");
    fflush(logfile);
  }

  //Copy everything into the right formats
  struct datapoint * data = (struct datapoint *) malloc( N * sizeof (struct datapoint) );
  struct gaussian * gaussians = (struct gaussian *) malloc (K * sizeof (struct gaussian) );

  bool noproj= (bool) noprojection;
  bool noweight= (bool) noweights;
  bool diagerrs= (bool) diagerrors;
  int ii, jj,dd1,dd2;
  for (ii = 0; ii != N; ++ii){
    data->ww = gsl_vector_alloc(dy);
    if ( ! noweight ) data->logweight = *(logweights++);
    if ( diagerrs ) data->SS = gsl_matrix_alloc(dy,1);
    else data->SS = gsl_matrix_alloc(dy,dy);
    if ( ! noproj ) data->RR = gsl_matrix_alloc(dy,d);
    for (dd1 = 0; dd1 != dy;++dd1)
      gsl_vector_set(data->ww,dd1,*(ydata++));
    if ( diagerrs)
      for (dd1 = 0; dd1 != dy; ++dd1)
	  gsl_matrix_set(data->SS,dd1,0,*(ycovar++));
    else
      for (dd1 = 0; dd1 != dy; ++dd1)
	for (dd2 = 0; dd2 != dy; ++dd2)
	  gsl_matrix_set(data->SS,dd1,dd2,*(ycovar++));
    if ( ! noproj )
      for (dd1 = 0; dd1 != dy; ++dd1)
	for (dd2 = 0; dd2 != d; ++dd2)
	  gsl_matrix_set(data->RR,dd1,dd2,*(projection++));
    else data->RR= NULL;
    ++data;
  }
  data -= N;
  ydata -= N*dy;
  if ( diagerrs ) ycovar -= N*dy;
  else ycovar -= N*dy*dy;
  if ( ! noproj ) projection -= N*dy*d;

  for (jj = 0; jj != K; ++jj){
    gaussians->mm = gsl_vector_alloc(d);
    gaussians->VV = gsl_matrix_alloc(d,d);
    gaussians->alpha = *(amp++);
    for (dd1 = 0; dd1 != d; ++dd1)
      gsl_vector_set(gaussians->mm,dd1,*(xmean++));
    for (dd1 = 0; dd1 != d; ++dd1)
      for (dd2 = 0; dd2 != d; ++dd2)
	gsl_matrix_set(gaussians->VV,dd1,dd2,*(xcovar++));
    ++gaussians;
  }
  gaussians -= K;
  amp -= K;
  xmean -= K*d;
  xcovar -= K*d*d;

  //Print the initial model parameters to the logfile
  int kk;
  if (keeplog){
    fprintf(logfile,"#\n#Using %i Gaussians and w = %f\n\n",K,w);
    fprintf(logfile,"#\n#Initial model parameters used:\n\n");
    for (kk=0; kk != K; ++kk){
      fprintf(logfile,"#Gaussian ");
      fprintf(logfile,"%i",kk);
      fprintf(logfile,"\n");
      fprintf(logfile,"#amp\t=\t");
      fprintf(logfile,"%f",(*gaussians).alpha);
      fprintf(logfile,"\n");
      fprintf(logfile,"#mean\t=\t");
      for (dd1=0; dd1 != d; ++dd1){
	fprintf(logfile,"%f",gsl_vector_get(gaussians->mm,dd1));
	if (dd1 < d-1) fprintf(logfile,"\t");
      }
      fprintf(logfile,"\n");
      fprintf(logfile,"#covar\t=\t");
      for (dd1=0; dd1 != d; ++dd1)
	fprintf(logfile,"%f\t",gsl_matrix_get(gaussians->VV,dd1,dd1));
      for (dd1=0; dd1 != d-1; ++dd1)
	for (dd2=dd1+1; dd2 != d; ++dd2){
	  fprintf(logfile,"%f\t",gsl_matrix_get(gaussians->VV,dd1,dd2));
	}
      ++gaussians;
      fprintf(logfile,"\n#\n");
    }
    gaussians -= K;
    fflush(logfile);
  }

  //Then run projected_gauss_mixtures
  proj_gauss_mixtures(data,N,gaussians,K,(bool *) fixamp,
		      (bool *) fixmean, (bool *) fixcovar,avgloglikedata,
		      tol,(long long int) maxiter, (bool) likeonly, w,
		      splitnmerge,keeplog,logfile,convlogfile,noproj,diagerrs,
		      noweight);

  //Print the final model parameters to the logfile
  if (keeplog){
    fprintf(logfile,"\n#Final model parameters obtained:\n\n");
    for (kk=0; kk != K; ++kk){
      fprintf(logfile,"#Gaussian ");
      fprintf(logfile,"%i",kk);
      fprintf(logfile,"\n");
      fprintf(logfile,"#amp\t=\t");
      fprintf(logfile,"%f",(*gaussians).alpha);
      fprintf(logfile,"\n");
      fprintf(logfile,"#mean\t=\t");
      for (dd1=0; dd1 != d; ++dd1){
	fprintf(logfile,"%f",gsl_vector_get(gaussians->mm,dd1));
	if (dd1 < d-1) fprintf(logfile,"\t");
      }
      fprintf(logfile,"\n");
      fprintf(logfile,"#covar\t=\t");
      for (dd1=0; dd1 != d; ++dd1)
	fprintf(logfile,"%f\t",gsl_matrix_get(gaussians->VV,dd1,dd1));
      for (dd1=0; dd1 != d-1; ++dd1)
	for (dd2=dd1+1; dd2 != d; ++dd2){
	  fprintf(logfile,"%f\t",gsl_matrix_get(gaussians->VV,dd1,dd2));
	}
      ++gaussians;
      fprintf(logfile,"\n#\n");
    }
    gaussians -= K;
    fflush(logfile);
  }

  //Then update the arrays given to us by IDL
  for (jj = 0; jj != K; ++jj){
    *(amp++) = gaussians->alpha;
    for (dd1 = 0; dd1 != d; ++dd1)
      *(xmean++) = gsl_vector_get(gaussians->mm,dd1);
    for (dd1 = 0; dd1 != d; ++dd1)
      for (dd2 = 0; dd2 != d; ++dd2)
	*(xcovar++) = gsl_matrix_get(gaussians->VV,dd1,dd2);
    ++gaussians;
  }
  gaussians -= K;
  amp -= K;
  xmean -= K*d;
  xcovar -= K*d*d;

  //And free any memory we allocated
  for (ii = 0; ii != N; ++ii){
    gsl_vector_free(data->ww);
    gsl_matrix_free(data->SS);
    if ( ! noproj )  gsl_matrix_free(data->RR);
    ++data;
  }
  data -= N;
  free(data);

  for (jj = 0; jj != K; ++jj){
    gsl_vector_free(gaussians->mm);
    gsl_matrix_free(gaussians->VV);
    ++gaussians;
  }
  gaussians -= K;
  free(gaussians);

  if (keeplog){
    fclose(logfile);
    fclose(convlogfile);
  }

  return 0;
}


text:
/**
 * @copyright (c) 2017 King Abdullah University of Science and Technology (KAUST).
 *                     All rights reserved.
 **/
/**
 * @file codelet_zhagdm.c
 *
 * Codelet for generating dense matrix from a problem determined according to current global setting of HiCMA library.
 *
 * HiCMA is a software package provided by King Abdullah University of Science and Technology (KAUST)
 *
 * @version 0.1.1
 * @author Kadir Akbudak
 * @date 2018-11-08
 **/
#include "morse.h"
#include "runtime/starpu/chameleon_starpu.h"

#include "runtime/starpu/runtime_codelets.h"
ZCODELETS_HEADER(hagdm)

#include <assert.h>
#include <stdio.h>
#include <sys/time.h>//FIXME for gettimeofday

#include "hicma.h"
#include "starsh.h"
#include "starsh-spatial.h"
#include "starsh-randtlr.h"
#ifdef MKL
  #include <mkl.h>
  #include <mkl_lapack.h>
  //#pragma message("MKL is used")
#else
  #ifdef ARMPL
    #include <armpl.h>
  #else
    #include <cblas.h>
  #endif
  #ifdef LAPACKE_UTILS
    #include <lapacke_utils.h>
  #endif
  #include <lapacke.h>
  //#pragma message("MKL is NOT used")
#endif

//#warning "An experimental feature is enabled!!!"
extern int steal_lrtile;

extern void _printmat(double * A, int m, int n, int ld);

void zhagdm(
        int nrows_Dense,
        int ncols_Dense,
        double *Dense,
        int ld_Dense,
        int tile_row_index,
        int tile_col_index,
        int A_mt
        )
{
    if(steal_lrtile == 1 && tile_row_index == tile_col_index){
        if(tile_row_index == A_mt-1) { // steal tile above
            tile_row_index -= 1;
        } else { // still tile below
            tile_row_index += 1;
        }
    }
    struct timeval tvalBefore, tvalAfter;
    gettimeofday (&tvalBefore, NULL);
    STARSH_cluster *RC = HICMA_get_starsh_format()->row_cluster, *CC = RC;
    void *RD = RC->data, *CD = RD;
    HICMA_get_starsh_format()->problem->kernel(nrows_Dense, ncols_Dense,
            RC->pivot+RC->start[tile_row_index],
            CC->pivot+CC->start[tile_col_index],
            RD, CD, Dense, ld_Dense);
}
/**
 * HICMA_TASK_zhagdm - Generate dense matrix from a problem determined according to current global setting of HiCMA library
 */
void HICMA_TASK_zhagdm( const MORSE_option_t *options,
                        int nrows_Dense, int ncols_Dense,
                        const MORSE_desc_t *Dense,
                        int ld_Dense,
                        int tile_row_index,
                        int tile_col_index,
                        int A_mt
                        )
{
    struct starpu_codelet *codelet = &cl_zhagdm;
    void (*callback)(void*) = NULL;

    MORSE_BEGIN_ACCESS_DECLARATION;
    MORSE_ACCESS_W(Dense, tile_row_index, tile_col_index);
    MORSE_END_ACCESS_DECLARATION;
    starpu_insert_task(
            starpu_mpi_codelet(codelet),
            STARPU_VALUE,    &nrows_Dense,                   sizeof(int),
            STARPU_VALUE,    &ncols_Dense,                   sizeof(int),
            STARPU_W,         RTBLKADDR(Dense, double, tile_row_index, tile_col_index),
            STARPU_VALUE,  &ld_Dense,                        sizeof(int),
            STARPU_VALUE,   &tile_row_index,                 sizeof(int),
            STARPU_VALUE,   &tile_col_index,                 sizeof(int),
            STARPU_VALUE,   &A_mt,                           sizeof(int),
            STARPU_PRIORITY,    options->priority,
            STARPU_CALLBACK,    callback,
#if defined(CHAMELEON_CODELETS_HAVE_NAME)
            STARPU_NAME, "zhagdm",
#endif
            0);
}
/**
 * cl_zhagdm_cpu_func - Generate a tile for random matrix.
 */

#if !defined(CHAMELEON_SIMULATION)
static void cl_zhagdm_cpu_func(void *descr[], void *cl_arg)
{
    int nrows_Dense;
    int ncols_Dense;
    int ld_Dense;
    int tile_row_index;
    int tile_col_index;
    int A_mt;
    int maxrank;
    double *Dense;

    Dense = (double *)STARPU_MATRIX_GET_PTR(descr[0]);
    starpu_codelet_unpack_args(cl_arg, &nrows_Dense, &ncols_Dense, &ld_Dense, &tile_row_index, &tile_col_index, &A_mt);
    zhagdm(
            nrows_Dense,
            ncols_Dense,
            Dense,
            ld_Dense,
            tile_row_index,
            tile_col_index,
            A_mt
            );
}
#endif /* !defined(CHAMELEON_SIMULATION) */

/*
 * Codelet definition
 */
CODELETS_CPU(zhagdm, 1, cl_zhagdm_cpu_func)

ZCODELETS_HEADER(hagdmi)
/**
 * HICMA_TASK_zhagdmi - Generate dense matrix from a problem determined according to current global setting of HiCMA library
 * This function takes indices of tiles of problem.
 */
void HICMA_TASK_zhagdmi( const MORSE_option_t *options,
                        int nrows_Dense, int ncols_Dense,
                        const MORSE_desc_t *Dense,
                        int ld_Dense,
                        int tile_row_index,
                        int tile_col_index,
                        int problem_row_index,
                        int problem_col_index
                        )
{
    struct starpu_codelet *codelet = &cl_zhagdmi;
    void (*callback)(void*) = NULL;

    MORSE_BEGIN_ACCESS_DECLARATION;
    MORSE_ACCESS_W(Dense, tile_row_index, tile_col_index);
    MORSE_END_ACCESS_DECLARATION;
    starpu_insert_task(
            starpu_mpi_codelet(codelet),
            STARPU_VALUE,    &nrows_Dense,                   sizeof(int),
            STARPU_VALUE,    &ncols_Dense,                   sizeof(int),
            STARPU_W,         RTBLKADDR(Dense, double, tile_row_index, tile_col_index),
            STARPU_VALUE,  &ld_Dense,                        sizeof(int),
            STARPU_VALUE,   &tile_row_index,                 sizeof(int),
            STARPU_VALUE,   &tile_col_index,                 sizeof(int),
            STARPU_VALUE,   &problem_row_index,                 sizeof(int),
            STARPU_VALUE,   &problem_col_index,                 sizeof(int),
            STARPU_PRIORITY,    options->priority,
            STARPU_CALLBACK,    callback,
#if defined(CHAMELEON_CODELETS_HAVE_NAME)
            STARPU_NAME, "zhagdm",
#endif
            0);
}

/** cl_zhagdm_cpu_func - Generate a tile for random matrix.
 * This function takes indices of tiles of problem.
 */

#if !defined(CHAMELEON_SIMULATION)
static void cl_zhagdmi_cpu_func(void *descr[], void *cl_arg)
{
    int nrows_Dense;
    int ncols_Dense;
    int ld_Dense;
    int tile_row_index;
    int tile_col_index;
    int maxrank;
    double *Dense;
    int problem_row_index;
    int problem_col_index;

    Dense = (double *)STARPU_MATRIX_GET_PTR(descr[0]);
    starpu_codelet_unpack_args(cl_arg, &nrows_Dense, &ncols_Dense, &ld_Dense, &tile_row_index, &tile_col_index, &problem_row_index, &problem_col_index);
    zhagdm(
            nrows_Dense,
            ncols_Dense,
            Dense,
            ld_Dense,
            problem_row_index,
            problem_col_index, -1
            );
}
#endif /* !defined(CHAMELEON_SIMULATION) */

/*
 * Codelet definition
 */
CODELETS_CPU(zhagdmi, 1, cl_zhagdmi_cpu_func)