text:
#ifndef LADMM_H
#define LADMM_H

#include "Matrix.h"

#include <string>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <stdio.h>      /* printf */
#include <time.h>
#include <fstream>
#include <algorithm>
#include <iomanip>
#include <ctime>
#include <sstream>
//#include "cmd_line.h"

//This class implements the method LADMM

/*
The optimization problem to solve is:

min \sum_i f^i(A_ix)+ \sum_i h^i(y^i)+ g(x)
s.t. Mx= y
Assumption 1: For each i, f_i is smooth, g(x) is seperable
*/

template<typename L, typename D>
class LADMM
{
private:

  std::vector<D> Ax;
  std::vector<D> old_Ax;

  std::vector<D> Mx;
  std::vector<D> old_Mx;

  std::vector<D> Mty;
  std::vector<D> old_Mty;

  std::vector<D> Mtlambda;
  std::vector<D> old_Mtlambda;

  std::vector<D> gradient;

  std::vector<D> MtMx;

protected:
  Matrix<L, D> data_A;

  Matrix<L, D> data_M;

  std::vector<D> x;

  std::vector<D> old_x;

  std::vector<D> y;

  std::vector<D> old_y;

  std::vector<D> lambda;

  std::vector<D> old_lambda;

  D beta;

  L m_1;

  L m_2;

  L m_3;

  D rho;

  D tau;

  D sigma;

  D L_phi;

  D function_value;

  D infeas;

  L print_every_N_ADMM;

  D running_time_ADMM;

  L nb_outer_iters;

  ofstream samp_ADMM;

public:
  D lambda_f;

  D mu_g;

  D lambda1;

  D lambda2;

  D L_h;

  virtual inline D value_of_f_j(D, L){return D(NULL);}
  virtual inline D value_of_h_j(D, L){return D(NULL);}
  virtual inline D gradient_of_f_j(D, L){return D(NULL);}
  virtual inline D prox_of_h_j(D,D, L){return D(NULL);}
  virtual inline D value_of_g_j(D, L){return D(NULL);}
  virtual inline D prox_of_g_j(D, D, L){return D(NULL);}
  virtual inline void set_matrix_M(){}
  virtual inline void set_matrix_A(){}

/*
  LADMM(const char* matrix_file, const char* matrix_file2)
  : Primal_Dual_LOOPLESS_Katyusha0<L,D>(),data_A(matrix_file), data_M(matrix_file2)
  {
  	this->matrix_merge(data_A,data_M);
    this->gamma=1;
  }
 */
  inline void set_L_phi(){
  	if (data_A.nsamples== 0){
  		L_phi= 0;
	  }
	else{
  		L_phi= compute_lambda_max_A(10);
  	}
  }

  D compute_lambda_max_A(L K){

    std::vector<D> bk(data_A.nfeatures);
    for (L j=0;j<data_A.nfeatures;j++)
    {
      bk[j]=1;
    }
    std::vector<D> yk(data_A.nsamples);
    D normk;
    D tmp;
    for(L kk=0;kk<K;kk++){
      for (L i=0;i<data_A.nsamples;i++){
        tmp=0;
        for (L k = data_A.ptr[i]; k < data_A.ptr[i + 1];	k++)
        {
          L j=data_A.row_idx[k];
          tmp+=data_A.A[k]*bk[j];
        }
        yk[i]=tmp;
      }
      normk=0;
      for (L j=0;j<data_A.nfeatures;j++){
        bk[j]=0;
        for (L k = data_A.ptr_t[j]; k < data_A.ptr_t[j + 1];	k++)
        {
          L i=data_A.col_idx[k];
          bk[j]+=data_A.A_t[k]*yk[i]*lambda_f;
        }
        normk+=bk[j]*bk[j];
      }
      normk=sqrt(normk);
      for (L j=0;j<data_A.nfeatures;j++)
      {bk[j]=bk[j]/normk; }
    }
    cout<<endl;
    D res=0;
    normk=0;
    for (L i=0;i<data_A.nsamples;i++){
      tmp=0;
      for (L k = data_A.ptr[i]; k < data_A.ptr[i + 1];	k++)
      {
        L j=data_A.row_idx[k];
        tmp+=data_A.A[k]*bk[j];
      }
      yk[i]=tmp;
      normk+=yk[i]*yk[i];
    }
    std::vector<D> bk2(data_A.nfeatures);
    for (L j=0;j<data_A.nfeatures;j++){
      bk2[j]=0;
      for (L k = data_A.ptr_t[j]; k < data_A.ptr_t[j + 1];	k++)
      {
        L i=data_A.col_idx[k];
        bk2[j]+=data_A.A_t[k]*yk[i]*lambda_f;
      }
    }
    for (L j=0;j<data_A.nfeatures;j++)
    res+=bk2[j]*bk[j];
    return res;
  }

D compute_lambda_max_M(L K){

    std::vector<D> bk(data_M.nfeatures);
    for (L j=0;j<data_M.nfeatures- 1;j++)
    {
      bk[j]=1;
    }
    bk[data_M.nfeatures- 1]= 2;
    std::vector<D> yk(data_M.nsamples);
    D normk;
    D tmp;
    for(L kk=0;kk<K;kk++){
      for (L i=0;i<data_M.nsamples;i++){
        tmp=0;
        for (L k = data_M.ptr[i]; k < data_M.ptr[i + 1];	k++)
        {
          L j=data_M.row_idx[k];
          tmp+=data_M.A[k]*bk[j];
        }
        yk[i]=tmp;
      }
      normk=0;
      for (L j=0;j<data_M.nfeatures;j++){
        bk[j]=0;
        for (L k = data_M.ptr_t[j]; k < data_M.ptr_t[j + 1];	k++)
        {
          L i=data_M.col_idx[k];
          bk[j]+=data_M.A_t[k]*yk[i];
        }
        normk+=bk[j]*bk[j];
      }
      normk=sqrt(normk);
      for (L j=0;j<data_M.nfeatures;j++)
      {bk[j]=bk[j]/normk; }
    }
    D res=0;
    normk=0;
    for (L i=0;i<data_M.nsamples;i++){
      tmp=0;
      for (L k = data_M.ptr[i]; k < data_M.ptr[i + 1];	k++)
      {
        L j=data_M.row_idx[k];
        tmp+=data_M.A[k]*bk[j];
      }
      yk[i]=tmp;
      normk+=yk[i]*yk[i];
    }
    std::vector<D> bk2(data_M.nfeatures);
    for (L j=0;j<data_M.nfeatures;j++){
      bk2[j]=0;
      for (L k = data_M.ptr_t[j]; k < data_M.ptr_t[j + 1];	k++)
      {
        L i=data_M.col_idx[k];
        bk2[j]+=data_M.A_t[k]*yk[i];
      }
    }
    for (L j=0;j<data_M.nfeatures;j++)
    res+=bk2[j]*bk[j];
    return res;
  }

  L get_nb_features(){
    return data_A.get_d();
  }

 /* inline D get_lambda1(){return lambda1;}
  inline D get_lambda2(){return lambda2;}
*/

  void update_x(){
    for (L i= 0; i< m_1; i++){
    	gradient[i]= 0;
    	for (L k = data_A.ptr_t[i]; k < data_A.ptr_t[i + 1];k++)
    	{
      		L j=data_A.col_idx[k];
      		gradient[i]+= data_A.A_t[k]*gradient_of_f_j(Ax[j],j);
    	}
    	MtMx[i]= 0;
    	for (L k = data_M.ptr_t[i]; k < data_M.ptr_t[i + 1];k++)
    	{
      		L j=data_M.col_idx[k];
      		MtMx[i]+= data_M.A_t[k]*Mx[j];
    	}
    	x[i]= prox_of_g_j(x[i]- (Mtlambda[i]+ gradient[i]+ beta*MtMx[i]- beta*Mty[i])/(beta*tau+ L_phi), beta*tau+ L_phi, i);
	}
	for (L i= 0; i< m_2; i++){
		Ax[i]= compute_AiTx(i);
	}
	for (L i= 0; i< m_3; i++){
		Mx[i]= compute_MiTx(i);
	}
  }

  void update_y(){
  	for (L i=0; i< m_3; i++){
  		y[i]= prox_of_h_j(y[i]+ (lambda[i]+ beta*Mx[i]- beta*y[i])/(beta*sigma), beta*sigma, i);
	  }
	for (L i=0; i< m_1; i++){
		Mty[i]= compute_MTiTy(i);
	}
  }

  D compute_AiTx(L i){
    D res=0;
    for (L k = data_A.ptr[i]; k < data_A.ptr[i + 1];k++)
    {
      L j=data_A.row_idx[k];
      res+= data_A.A[k]*x[j];
    }
    return res;
  }

  D compute_MiTx(L i){
    D res=0;
    for (L k = data_M.ptr[i]; k < data_M.ptr[i + 1];k++)
    {
      L j=data_M.row_idx[k];
      res+= data_M.A[k]*x[j];
    }
    return res;
  }

  D compute_MTiTy(L i){
    D res=0;
    for (L k = data_M.ptr_t[i]; k < data_M.ptr_t[i + 1];k++)
    {
      L j=data_M.col_idx[k];
      res+= data_M.A_t[k]*y[j];
    }
    return res;
  }

  D compute_MTiTlambda(L i){
    D res=0;
    for (L k = data_M.ptr_t[i]; k < data_M.ptr_t[i + 1];k++)
    {
      L j=data_M.col_idx[k];
      res+= data_M.A_t[k]*lambda[j];
    }
    return res;
  }

  void compute_function_value(){
    D res= 0;
    for (L i= 0; i< data_A.nfeatures; i++){
      res+= value_of_g_j(x[i],i);
    }
    for (L i= 0; i< data_M.nsamples; i++){
      res+= value_of_h_j(Mx[i],i);
    }
    for (L i= 0; i<data_A.nsamples; i++){
      res+= value_of_f_j(Ax[i],i);
    }
    function_value= res;
  }

  void compute_infeasibility(){
  	D res = 0;
  	for (L i= 0; i< data_M.nsamples; i++){
      res+= (Mx[i]- data_M.b[i])*(Mx[i]- data_M.b[i]);
    }
    infeas= sqrt(res);
  }

void update_lambda(){
  for(L j=0;j<m_3;j++)
  {
    lambda[j]= lambda[j]+ beta*(Mx[j]- y[j]);
  }
  for (L i=0; i< m_1; i++){
		Mtlambda[i]= compute_MTiTlambda(i);
	}
}

void Initialize(D beta_0, D val_rho, vector<D> & x0,vector<D> & y0, vector<D> & lambda0){
  cout<<"start initializing"<<endl;
  set_matrix_M();
  set_matrix_A();

  m_1=data_A.get_d();
  m_2=data_A.get_n();
  m_3=data_M.get_n();
  cout<<"m_1="<<m_1<<endl;
  cout<<"m_2="<<m_2<<endl;
  cout<<"m_3="<<m_3<<endl;

  beta=beta_0;
  //tau= data_A.nsamples;

  tau= 1.02*compute_lambda_max_M(10);
  sigma= 1;

  rho=val_rho;

  x.resize(m_1,0);
  old_x.resize(m_1,0);
  y.resize(m_3,0);
  old_y.resize(m_3,0);
  lambda.resize(m_3,0);
  old_lambda.resize(m_3,0);
  for(L i=0;i<m_1;i++){
    x[i]=x0[i];
  }
  for(L j=0;j<m_3;j++){
    y[j]=y0[j];
  }
  for(L j=0;j<m_3;j++){
    lambda[j]=lambda0[j];
  }

  Ax.clear();
  Ax.resize(m_2,0);
  Mx.clear();
  Mx.resize(m_3,0);
  Mty.clear();
  Mty.resize(m_1,0);
  Mtlambda.clear();
  Mtlambda.resize(m_1,0);
  gradient.clear();
  gradient.resize(m_1,0);
  MtMx.clear();
  MtMx.resize(m_1,0);
  L_phi= 0;
  set_L_phi();
  cout<< "L_phi= "<< L_phi << " beta= "<< beta<< " rho= "<< rho<< " tau= "<< tau<< " sigma= "<< sigma<<endl;
}

void reset_everything(){
  beta*=rho;

}

inline void compute_and_record_res(){
  if(nb_outer_iters%print_every_N_ADMM==0){
    compute_function_value();
    compute_infeasibility();
    cout<<setprecision(9)<<"Iteration: "<<nb_outer_iters<<"; time="<<running_time_ADMM<< "; function value="<<function_value<< "; infeasibility="<< infeas<< endl;
    samp_ADMM<<setprecision(9)<<nb_outer_iters<<" "<<running_time_ADMM<<" "<< function_value<<" "<< infeas<< endl;
  }
}

void ADMM_solve_with_Linear(D beta_0, D val_rho,vector<D> & x0,vector<D> & y0, vector<D> & lambda0, L max_nb_outer, L p_N_1, string filename1, D time){
  Initialize(beta_0,val_rho, x0, y0, lambda0);
  nb_outer_iters=0;
  //string sampname2= ALGparam.data_dir +"/results/L_Katyusha_"+filename2;
  //filename1= ALGparam.data_dir +"/results/ADMM_"+filename1;
  filename1= "results/ADMM_"+filename1;
  samp_ADMM.open(filename1.c_str());
  running_time_ADMM=0;
  print_every_N_ADMM=p_N_1;
  compute_and_record_res();
  D start;
  D res_x, res_y, res_l;
  /*
  for(L i=0;i<m_3;i++){
    old_lambda[i]=lambda[i];
  }
  */
  while(nb_outer_iters<max_nb_outer){
    //rescale();
    for(L i=0;i<m_1;i++){
      old_x[i]=x[i];
    }
    for(L i=0;i<m_3;i++){
      old_y[i]=y[i];
    }
    for(L i=0;i<m_3;i++){
      old_lambda[i]=lambda[i];
    }
    start = std::clock();
    update_x();
    update_y();
    update_lambda();
    nb_outer_iters++;
    running_time_ADMM+=( std::clock() - start ) / (double) CLOCKS_PER_SEC;
    compute_and_record_res();
    /*
    res_x= 0;
    for(L i=0;i<m_1;i++){
      res_x+= (old_x[i]- x[i])*(old_x[i]- x[i]);
    }
    res_y= 0;
    for(L i=0;i<m_3;i++){
      res_y= (old_y[i]- y[i])*(old_y[i]- y[i]);
    }
    res_l= 0;
    for(L i=0;i<m_3;i++){
      res_l= (old_lambda[i]- lambda[i])*(old_lambda[i]- lambda[i]);
    }
    cout<< "res_x= "<< res_x<< " res_y= "<< res_y<< " res_l= "<< res_l<< endl;
    system("pause");
    */
    start = std::clock();
    reset_everything();
    running_time_ADMM+=( std::clock() - start ) / (double) CLOCKS_PER_SEC;
    if (running_time_ADMM> time){
      break;
    }
  }

}
};

#endif /* MIN_SMOOTH_CONVEX_H */


text:
/**
 *
 * @file qwrapper_spotrf.c
 *
 *  PLASMA core_blas quark wrapper
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Hatem Ltaief
 * @author Mathieu Faverge
 * @author Jakub Kurzak
 * @date 2010-11-15
 * @generated s Tue Jan  7 11:44:56 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 **/
void QUARK_CORE_spotrf(Quark *quark, Quark_Task_Flags *task_flags,
                       PLASMA_enum uplo, int n, int nb,
                       float *A, int lda,
                       PLASMA_sequence *sequence, PLASMA_request *request,
                       int iinfo)
{
    DAG_CORE_POTRF;
    QUARK_Insert_Task(quark, CORE_spotrf_quark, task_flags,
        sizeof(PLASMA_enum),                &uplo,      VALUE,
        sizeof(int),                        &n,         VALUE,
        sizeof(float)*nb*nb,    A,                 INOUT,
        sizeof(int),                        &lda,       VALUE,
        sizeof(PLASMA_sequence*),           &sequence,  VALUE,
        sizeof(PLASMA_request*),            &request,   VALUE,
        sizeof(int),                        &iinfo,     VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_spotrf_quark = PCORE_spotrf_quark
#define CORE_spotrf_quark PCORE_spotrf_quark
#endif
void CORE_spotrf_quark(Quark *quark)
{
    PLASMA_enum uplo;
    int n;
    float *A;
    int lda;
    PLASMA_sequence *sequence;
    PLASMA_request *request;
    int iinfo;

    int info;

    quark_unpack_args_7(quark, uplo, n, A, lda, sequence, request, iinfo);
    info = LAPACKE_spotrf_work(
        LAPACK_COL_MAJOR,
        lapack_const(uplo),
        n, A, lda);
    if (sequence->status == PLASMA_SUCCESS && info != 0)
        plasma_sequence_flush(quark, sequence, request, iinfo+info);
}

ext: c


text:
/*
Copyright (C) 2019-2020 JingWeiZhangHuai <jingweizhanghuai@163.com>
Licensed under the Apache License, Version 2.0; 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.
*/

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

#include <cblas.h>
#include "morn_tensor.h"

void DeconvTensorToMatData(MTensor *tns,int bc,float *mdata,int knl_height,int knl_width,int y_stride,int x_stride)
{
    int height = tns->height;
    int width  = tns->width;
    int channel= tns->channel;

    int out_width = width*x_stride;
    int out_height=height*y_stride;
    int mwidth = knl_height*knl_width*channel+1;
    int mheight= out_height*out_width;

    memset(mdata,0,mwidth*mheight*sizeof(float));

    float *tdata = tns->data[bc];
    int tsize = tns->height*tns->width;

    int i,j,c;
    for(j=0;j<mheight;j++)
    {
        int n=j/out_width-knl_height/2;
        int m=j%out_width-knl_width/2;
        for(i=0;i<mwidth-1;i+=channel)
        {
            int h= i/(knl_width*channel)         +n;if(h%y_stride!=0) continue;h=h/y_stride;
            int w=(i%(knl_width*channel))/channel+m;if(w%x_stride!=0) continue;w=w/x_stride;

            if(h<0)h=0;else if(h>=height)h=height-1;
            if(w<0)w=0;else if(w>= width)w= width-1;

            for(c=0;c<channel;c++)
                mdata[(j*mwidth)+i+c]=tdata[c*tsize+h*width+w];
        }
        mdata[(j*mwidth)+mwidth-1]=1.0f;
    }
}

void DeconvMatDataToTensor(float *mdata,MTensor *tns,int bc,int knl_height,int knl_width,int y_stride,int x_stride)
{
    int height = tns->height;
    int width  = tns->width;
    int channel= tns->channel;

    int out_width = width*x_stride;
    int out_height=height*y_stride;

    int mwidth = knl_height*knl_width*channel+1;
    int mheight= out_height*out_width;

    float *tdata = tns->data[bc];
    int tsize = tns->height*tns->width;

    int i,j,c;
    for(j=0;j<mheight;j++)
    {
        int n=j/out_width-knl_height/2;
        int m=j%out_width-knl_width/2;
        for(i=0;i<mwidth-1;i+=channel)
        {
            int h= i/(knl_width*channel)         +n;if(h%y_stride!=0) continue;h=h/y_stride;
            int w=(i%(knl_width*channel))/channel+m;if(w%x_stride!=0) continue;w=w/x_stride;

            if(h<0)h=0;else if(h>=height)h=height-1;
            if(w<0)w=0;else if(w>= width)w= width-1;

            for(c=0;c<channel;c++)
                tdata[c*tsize+h*width+w]+=mdata[(j*mwidth)+i+c];
        }
        mdata[(j*mwidth)+mwidth-1]=1.0f;
    }
}

struct TensorDeconvPara
{
    MLayer *prev;

    int knl_num;
    int knl_height;
    int knl_width;

    int x_stride;
    int y_stride;

    int res_valid;

    float rate;
    float decay;
    float momentum;
};
void *mTensorDeconvPara(MFile *ini,char *name)
{
    struct TensorDeconvPara *para = (struct TensorDeconvPara *)mMalloc(sizeof(struct TensorDeconvPara));

    char *value = mINIRead(ini,name,"prev");
    para->prev = mNetworkLayer(ini,value);
    mException((para->prev == NULL),EXIT,"invalid prev");

    para->res_valid = (strcmp("Input",mLayerType(para->prev))!=0);

    value = mINIRead(ini,name,"knl_num");
    if(value != NULL) para->knl_num= atoi(value);else para->knl_num= 1;

    value = mINIRead(ini,name,"knl_height");
    if(value != NULL) para->knl_height= atoi(value);else para->knl_height= 1;

    value = mINIRead(ini,name,"knl_width");
    if(value != NULL) para->knl_width= atoi(value);else para->knl_width= 1;

    value = mINIRead(ini,name,"x_stride");
    if(value != NULL) para->x_stride= atoi(value);else para->x_stride= 1;

    value = mINIRead(ini,name,"y_stride");
    if(value != NULL) para->y_stride= atoi(value);else para->y_stride= 1;

    value = mINIRead(ini,name,"rate");
    if(value != NULL) para->rate = atof(value);
    else
    {
        value = mINIRead(ini,"para","rate");
        if(value != NULL) para->rate = atof(value);
        else              para->rate = 0.001;
    }

    value = mINIRead(ini,name,"decay");
    if(value != NULL) para->decay = atof(value);
    else
    {
        value = mINIRead(ini,"para","decay");
        if(value != NULL) para->decay = atof(value);
        else              para->decay = 0.01;
    }
    mException((para->decay<0.0f)||(para->decay>=1.0f),EXIT,"invalid para decay");

    value = mINIRead(ini,name,"momentum");
    if(value != NULL) para->momentum = atof(value);
    else
    {
        value = mINIRead(ini,"para","momentum");
        if(value != NULL) para->momentum = atof(value);
        else              para->momentum = 0.9;
    }

    return para;
}

struct HandleTensorDeconv
{
    float *mat;
    float *kernel;
    float *update;
};
void endTensorDeconv(void *info)
{
    struct HandleTensorDeconv *handle = (struct HandleTensorDeconv *)info;
    if(handle->mat   != NULL) mFree(handle->mat);
    if(handle->kernel!= NULL) mFree(handle->kernel);
    if(handle->update!= NULL) mFree(handle->update);
}
#define HASH_TensorDeconv 0x9087d39c
void TensorDeconvSet(MLayer *layer)
{
    if(layer->state != DFLT) return;
    struct TensorDeconvPara *para = (struct TensorDeconvPara *)(layer->para);
    MTensor *in = para->prev->tns;
    MTensor *res= para->prev->res;
    MTensor *out=layer->tns;

    MHandle *hdl=mHandle(out,TensorDeconv);
    struct HandleTensorDeconv *handle = (struct HandleTensorDeconv *)(hdl->handle);

    int out_height= in->height*para->y_stride;
    int out_width = in->width *para->x_stride;
    int mheight = (out_height*out_width);
    int mwidth = para->knl_height*para->knl_width*in->channel+1;
    int data_size = para->knl_num*mwidth;

    mTensorRedefine(out,in->batch,para->knl_num,out_height,out_width,NULL);
    if(morn_network_flag == MORN_TRAIN)
    {
        if(INVALID_TENSOR(res)) mTensorRedefine(res,in->batch,in->channel,in->height,in->width,in->data);
        else                    mTensorRedefine(res,in->batch,in->channel,in->height,in->width,NULL);

        if(handle->update != NULL) mFree(handle->update);
        handle->update =(float *)mMalloc(data_size*sizeof(float));
        memset(handle->update,0,data_size*sizeof(float));
    }

    if(handle->kernel !=NULL) mFree(handle->kernel);
    handle->kernel = (float *)mMalloc(data_size*sizeof(float));

    if(morn_network_parafile==NULL)
    {
        float scale = sqrt(2.0f/mwidth);
        for(int i=0;i<data_size;i++)
            handle->kernel[i] = scale*mNormalRand(0.0f,1.0f);
    }
    else
        mNetworkParaRead(layer,"kernel",handle->kernel,data_size*sizeof(float));

    if(handle->mat!=NULL) mFree(handle->mat);
    handle->mat = (float *)mMalloc(mheight*mwidth*sizeof(float));

    hdl->valid = 1;
}

void mTensorDeconvForward(MLayer *layer)
{
    mException(INVALID_POINTER(layer),EXIT,"invalid input");
    mException(strcmp("Deconv",mLayerType(layer)),EXIT,"invalid layer type");
    struct TensorDeconvPara *para = (struct TensorDeconvPara *)(layer->para);
    MTensor *in = para->prev->tns;
    MTensor *out=layer->tns;

    TensorDeconvSet(layer);

    MHandle *hdl=mHandle(out,TensorDeconv);
    struct HandleTensorDeconv *handle = (struct HandleTensorDeconv *)(hdl->handle);

    int mheight = (out->height*out->width);
    int mwidth = para->knl_height*para->knl_width*in->channel+1;

    float *kernel_data= handle->kernel;
    float *in_data = handle->mat;

    for(int b=0;b<in->batch;b++)
    {
        DeconvTensorToMatData(in,b,in_data,para->knl_height,para->knl_width,para->y_stride,para->x_stride);
        float *out_data = out->data[b];

        in_data[mwidth-1]=1.0f;

        cblas_sgemm(CblasRowMajor,CblasNoTrans,CblasTrans,
                    para->knl_num,mheight,mwidth,
                    1.0f,
                    kernel_data,mwidth,
                        in_data,mwidth,
               0.0f,   out_data,mheight);
    }

    layer->state = MORN_FORWARD;
}

void mTensorDeconvBackward(MLayer *layer)
{
    mException(INVALID_POINTER(layer),EXIT,"invalid input");
    mException(strcmp("Deconv",mLayerType(layer)),EXIT,"invalid layer type");
    struct TensorDeconvPara *para = (struct TensorDeconvPara *)(layer->para);

    MTensor *in = para->prev->tns;
    MTensor *res= para->prev->res;
    MTensor *out=layer->res;

    MHandle *hdl=mHandle(layer->tns,TensorDeconv);
    struct HandleTensorDeconv *handle = (struct HandleTensorDeconv *)(hdl->handle);
    mException((hdl->valid == 0),EXIT,"no forward operate");

    int mheight = (out->height*out->width);
    int mwidth = para->knl_height*para->knl_width*in->channel+1;

    float *kernel_data= handle->kernel;
    float *update_data= handle->update;
    float *    in_data= handle->mat;
    float *   res_data= handle->mat;

    mNetworkParaWrite(layer,"kernel",kernel_data,para->knl_num*mwidth*sizeof(float));

    for(int b=0;b<in->batch;b++)
    {
        DeconvTensorToMatData(in,b,in_data,para->knl_height,para->knl_width,para->y_stride,para->x_stride);
        float *out_data = out->data[b];

        in_data[mwidth-1]=1.0f;

        cblas_sgemm(CblasRowMajor,CblasNoTrans,CblasNoTrans,
                    para->knl_num,mwidth,mheight,
                    1.0f,
                       out_data,mheight,
                        in_data,mwidth,
                    (b==0)?para->momentum:1.0f,
                    update_data,mwidth);
    }
    cblas_saxpby(para->knl_num*mwidth,
                 (0.0f-(para->rate/(float)(in->batch))),update_data,1,
                 (1.0f-(para->decay*para->rate))       ,kernel_data,1);

    if(para->res_valid==0) return;

    if(para->prev->state == MORN_FORWARD)
    {
        for(int b=0;b<res->batch;b++)
            memset(res->data[b],0,in->height*in->width*in->channel*sizeof(float));
        para->prev->state = MORN_BACKWARD;
    }

    for(int b=0;b<in->batch;b++)
    {
        float *out_data = out->data[b];

        cblas_sgemm(CblasRowMajor,CblasTrans,CblasNoTrans,
                    mheight,mwidth,para->knl_num,
                    1.0f,
                       out_data,mheight,
                    kernel_data,mwidth,
                0.0,   res_data,mwidth);

        DeconvMatDataToTensor(res_data,res,b,para->knl_height,para->knl_width,para->y_stride,para->x_stride);
    }
}


text:
#include "spinodal_spect.h"
#include <gsl/gsl_rng.h>

// Generate array of random numbers between 0 and 1
double *rand_ZeroToOne(int Nx, int Ny, int seed, double *random_ZeroToOne_array) {
        const gsl_rng_type *T;
        gsl_rng *r;

        int i,
        NxNy = Nx * Ny;

        gsl_rng_env_setup();
        T = gsl_rng_default;
        r = gsl_rng_alloc(T);
        gsl_rng_set(r, seed);

        // Generate array of random numbers between 0 and 1
        for (i = 0; i < NxNy; i++) {
            random_ZeroToOne_array[i] = gsl_rng_uniform(r);
        }

        gsl_rng_free(r);

        return random_ZeroToOne_array;

    }


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

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_sf.h>
#include "test_sf.h"

#define PRINT(n) printf("%22.18g  %22.18g  %22.18g  %22.18g\n", F[n], Fp[n], G[n], Gp[n])

#define WKB_TOL (1.0e+04 * TEST_SQRT_TOL0)

int test_coulomb(void)
{
  gsl_sf_result r;
  int status = 0;
  int s = 0;

  char message_buff[2048];

  /* const int kmax = 20; */
  /* double F[kmax+1], Fp[kmax+1], G[kmax+1], Gp[kmax+1]; */
  gsl_sf_result F, Fp, G, Gp;
  double Fe, Ge;
  double lam_min;
  double lam_F;
  double eta, x;
  int k_G;

  TEST_SF(s, gsl_sf_hydrogenicR_1_e, (3.0, 2.0, &r),  0.025759948256148471036,  TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hydrogenicR_1_e, (3.0, 10.0, &r), 9.724727052062819704e-13, TEST_TOL1, GSL_SUCCESS);
  status += s;

  TEST_SF(s, gsl_sf_hydrogenicR_e, (4, 1, 3.0, 0.0, &r),  0.0,  TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hydrogenicR_e, (4, 0, 3.0, 2.0, &r), -0.03623182256981820062,  TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hydrogenicR_e, (4, 1, 3.0, 2.0, &r), -0.028065049083129581005, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hydrogenicR_e, (4, 2, 3.0, 2.0, &r),  0.14583027278668431009,  TEST_TOL0, GSL_SUCCESS);
  status += s;

  TEST_SF(s, gsl_sf_hydrogenicR_e, (100,  0, 3.0, 2.0, &r), -0.00007938950980052281367, TEST_TOL3, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hydrogenicR_e, (100, 10, 3.0, 2.0, &r),  7.112823375353605977e-12,  TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_hydrogenicR_e, (100, 90, 3.0, 2.0, &r),  5.845231751418131548e-245, TEST_TOL2, GSL_SUCCESS);
  status += s;

  lam_F = 0.0;
  k_G   = 0;
  eta = 1.0;
  x = 5.0;
  gsl_sf_coulomb_wave_FG_e(eta, x, lam_F, k_G, &F, &Fp, &G, &Gp, &Fe, &Ge);
  s = 0;
  message_buff[0] = 0;
  s += test_sf_check_result(message_buff,  F,  0.6849374120059439677, TEST_TOL3);
  s += test_sf_check_result(message_buff, Fp, -0.7236423862556063963, TEST_TOL3);
  s += test_sf_check_result(message_buff,  G, -0.8984143590920205487, TEST_TOL3);
  s += test_sf_check_result(message_buff, Gp, -0.5108047585190350106, TEST_TOL3);
  printf("%s", message_buff);
  gsl_test(s,"  gsl_sf_coulomb_wave_FG_e(1.0, 5.0, lam_F=0, lam_G=0)");
  status += s;

  lam_F = 10.0;
  k_G   = 2;
  eta = 1.0;
  x = 5.0;
  gsl_sf_coulomb_wave_FG_e(eta, x, lam_F, k_G, &F, &Fp, &G, &Gp, &Fe, &Ge);
  s = 0;
  message_buff[0] = 0;
  s += test_sf_check_result(message_buff,  F,  0.0006423773354915823698, TEST_TOL3);
  s += test_sf_check_result(message_buff, Fp,  0.0013299570958719702545, TEST_TOL3);
  s += test_sf_check_result(message_buff,  G,  33.27615734455096130,     TEST_TOL3);
  s += test_sf_check_result(message_buff, Gp, -45.49180102261540580,     TEST_TOL3);
  printf("%s", message_buff);
  gsl_test(s,"  gsl_sf_coulomb_wave_FG_e(1.0, 5.0, lam_F=10, lam_G=8)");
  status += s;

  lam_F = 4.0;
  k_G   = 2;
  eta = 50.0;
  x = 120.0;
  gsl_sf_coulomb_wave_FG_e(eta, x, lam_F, k_G, &F, &Fp, &G, &Gp, &Fe, &Ge);
  s = 0;
  message_buff[0] = 0;
  s += test_sf_check_result(message_buff,  F,  0.0735194711823798495, TEST_TOL3);
  s += test_sf_check_result(message_buff, Fp,  0.6368149124126783325, TEST_TOL3);
  /*
  s += test_sf_check_result(message_buff,  G,  , TEST_TOL5);
  s += test_sf_check_result(message_buff, Gp, , TEST_TOL5);
  */
  printf("%s", message_buff);
  gsl_test(s,"  gsl_sf_coulomb_wave_FG_e(50.0, 120.0, lam_F=4, lam_G=2)");
  status += s;

  lam_F = 0.0;
  k_G = 0;
  eta = -1000.0;
  x = 1.0;
  gsl_sf_coulomb_wave_FG_e(eta, x, lam_F, k_G, &F, &Fp, &G, &Gp, &Fe, &Ge);
  s = 0;
  message_buff[0] = 0;
  s += test_sf_check_result(message_buff,  F,  9.68222518991341e-02, TEST_TOL3);
  s += test_sf_check_result(message_buff, Fp,  5.12063396274631e+00, TEST_TOL3);
  s += test_sf_check_result(message_buff,  G,  1.13936784379472e-01, TEST_TOL3);
  s += test_sf_check_result(message_buff, Gp, -4.30243486522438e+00, TEST_TOL3);
  printf("%s", message_buff);
  gsl_test(s, "  gsl_sf_coulomb_wave_FG_e(-1000.0, 1.0, lam_F=0, lam_G=0)");
  status += s;

  lam_min = 0.0;
  eta = -50.0;
  x = 5.0;
  gsl_sf_coulomb_wave_FG_e(eta, x, lam_F, k_G, &F, &Fp, &G, &Gp, &Fe, &Ge);
  s = 0;
  message_buff[0] = 0;
  s += test_sf_check_result(message_buff,  F,  1.52236975714236e-01, TEST_TOL3);
  s += test_sf_check_result(message_buff, Fp,  2.03091041166137e+00, TEST_TOL3);
  s += test_sf_check_result(message_buff,  G,  4.41680690236251e-01, TEST_TOL3);
  s += test_sf_check_result(message_buff, Gp, -6.76485374766869e-01, TEST_TOL3);
  printf("%s", message_buff);
  gsl_test(s, "  gsl_sf_coulomb_wave_FG_e(-50.0, 5.0, lam_F=0, lam_G=0)");
  status += s;

  lam_min = 0.0;
  eta = -50.0;
  x = 1000.0;
  gsl_sf_coulomb_wave_FG_e(eta, x, lam_F, k_G, &F, &Fp, &G, &Gp, &Fe, &Ge);
  s = 0;
  message_buff[0] = 0;
  s += test_sf_check_result(message_buff,  F, -0.2267212182760888523, TEST_TOL3);
  s += test_sf_check_result(message_buff, Fp, -0.9961306810018401525, TEST_TOL3);
  s += test_sf_check_result(message_buff,  G, -0.9497684438900352186, TEST_TOL3);
  s += test_sf_check_result(message_buff, Gp,  0.2377656295411961399, TEST_TOL3);
  printf("%s", message_buff);
  gsl_test(s, "  gsl_sf_coulomb_wave_FG_e(-50.0, 1000.0, lam_F=0, lam_G=0)");
  status += s;

  lam_F = 10.0;
  k_G = 0;
  eta = -50.0;
  x = 5.0;
  gsl_sf_coulomb_wave_FG_e(eta, x, lam_F, k_G, &F, &Fp, &G, &Gp, &Fe, &Ge);
  s = 0;
  message_buff[0] = 0;
  s += test_sf_check_result(message_buff,  F, -3.681143602184922e-01, TEST_TOL3);
  s += test_sf_check_result(message_buff, Fp,  1.338467510317215e+00, TEST_TOL3);
  s += test_sf_check_result(message_buff,  G,  3.315883246109351e-01, TEST_TOL3);
  s += test_sf_check_result(message_buff, Gp,  1.510888628136180e+00, TEST_TOL3);
  printf("%s", message_buff);
  gsl_test(s, "  gsl_sf_coulomb_wave_FG_e(-50.0, 5.0, lam_F=10, lam_G=10)");
  status += s;

  lam_F = 0.0;
  k_G = 0;
  eta = -4.0;
  x = 5.0;
  gsl_sf_coulomb_wave_FG_e(eta, x, lam_F, k_G, &F, &Fp, &G, &Gp, &Fe, &Ge);
  s = 0;
  message_buff[0] = 0;
  s += test_sf_check_result(message_buff,  F,  4.078627230056172e-01, TEST_TOL3);
  s += test_sf_check_result(message_buff, Fp,  1.098212336357310e+00, TEST_TOL3);
  s += test_sf_check_result(message_buff,  G,  6.743270353832442e-01, TEST_TOL3);
  s += test_sf_check_result(message_buff, Gp, -6.361104272804447e-01, TEST_TOL3);
  printf("%s", message_buff);
  gsl_test(s, "  gsl_sf_coulomb_wave_FG_e(-4.0, 5.0, lam_F=0, lam_G=0");
  status += s;

  lam_F = 3.0;
  k_G = 0;
  eta = -4.0;
  x = 5.0;
  gsl_sf_coulomb_wave_FG_e(eta, x, lam_F, k_G, &F, &Fp, &G, &Gp, &Fe, &Ge);
  s = 0;
  message_buff[0] = 0;
  s += test_sf_check_result(message_buff,  F, -2.568630935581323e-01, TEST_TOL3);
  s += test_sf_check_result(message_buff, Fp,  1.143229422014827e+00, TEST_TOL3);
  s += test_sf_check_result(message_buff,  G,  7.879899223927996e-01, TEST_TOL3);
  s += test_sf_check_result(message_buff, Gp,  3.859905878106713e-01, TEST_TOL3);
  printf("%s", message_buff);
  gsl_test(s, "  gsl_sf_coulomb_wave_FG_e(-4.0, 5.0, lam_F=3, lam_G=3");
  status += s;

  lam_F = 0.0;
  k_G = 0;
  eta = 1.0;
  x = 2.0;
  gsl_sf_coulomb_wave_FG_e(eta, x, lam_F, k_G, &F, &Fp, &G, &Gp, &Fe, &Ge);
  s = 0;
  message_buff[0] = 0;
  s += test_sf_check_result(message_buff,  F,  6.61781613832681e-01, TEST_TOL3);
  s += test_sf_check_result(message_buff, Fp,  4.81557455709949e-01, TEST_TOL3);
  s += test_sf_check_result(message_buff,  G,  1.27577878476828e+00, TEST_TOL3);
  s += test_sf_check_result(message_buff, Gp, -5.82728813097184e-01, TEST_TOL3);
  printf("%s", message_buff);
  gsl_test(s, "  gsl_sf_coulomb_wave_FG_e(1.0, 2.0, lam_F=0, lam_G=0)");
  status += s;

  lam_F = 0.0;
  k_G = 0;
  eta = 1.0;
  x = 0.5;
  gsl_sf_coulomb_wave_FG_e(eta, x, lam_F, k_G, &F, &Fp, &G, &Gp, &Fe, &Ge);
  s = 0;
  message_buff[0] = 0;
  s += test_sf_check_result(message_buff,  F,  0.08315404535022023302, TEST_TOL3);
  s += test_sf_check_result(message_buff, Fp,  0.22693874616222787568, TEST_TOL3);
  s += test_sf_check_result(message_buff,  G,  3.1060069279548875140,  TEST_TOL3);
  s += test_sf_check_result(message_buff, Gp, -3.549156038719924236,   TEST_TOL3);
  printf("%s", message_buff);
  gsl_test(s, "  gsl_sf_coulomb_wave_FG_e(1.0, 0.5, lam_F=0, lam_G=0)");
  status += s;

  lam_F = 0.5;
  k_G = 0;
  eta = 1.0;
  x = 0.5;
  gsl_sf_coulomb_wave_FG_e(eta, x, lam_F, k_G, &F, &Fp, &G, &Gp, &Fe, &Ge);
  s = 0;
  message_buff[0] = 0;
  s += test_sf_check_result(message_buff,  F,  0.04049078073829290935, TEST_TOL3);
  s += test_sf_check_result(message_buff, Fp,  0.14194939168094778795, TEST_TOL3);
  s += test_sf_check_result(message_buff,  G,  4.720553853049677897,   TEST_TOL3);
  s += test_sf_check_result(message_buff, Gp, -8.148033852319180005,   TEST_TOL3);
  printf("%s", message_buff);
  gsl_test(s, "  gsl_sf_coulomb_wave_FG_e(1.0, 0.5, lam_F=0.5, lam_G=0.5)");
  status += s;

  lam_F = 0.1;
  k_G = 0;
  eta = 1.0;
  x = 0.5;
  gsl_sf_coulomb_wave_FG_e(eta, x, lam_F, k_G, &F, &Fp, &G, &Gp, &Fe, &Ge);
  s = 0;
  message_buff[0] = 0;
  s += test_sf_check_result(message_buff,  F,  0.07365466672379703418, TEST_TOL5);
  s += test_sf_check_result(message_buff, Fp,  0.21147121807178518647, TEST_TOL5);
  s += test_sf_check_result(message_buff,  G,  3.306705446241024890, TEST_TOL5);
  s += test_sf_check_result(message_buff, Gp, -4.082931670935696644, TEST_TOL5);
  printf("%s", message_buff);
  gsl_test(s, "  gsl_sf_coulomb_wave_FG_e(1.0, 0.5, lam_F=0.1, lam_G=0.1)");
  status += s;

  lam_F = 0.0;
  k_G = 0;
  eta = 8.0;
  x = 1.05;
  gsl_sf_coulomb_wave_FG_e(eta, x, lam_F, k_G, &F, &Fp, &G, &Gp, &Fe, &Ge);
  s = 0;
  message_buff[0] = 0;
  s += test_sf_check_result(message_buff,  F,  9.882706082810274357e-09, TEST_TOL5);
  s += test_sf_check_result(message_buff, Fp,  4.005167028235547770e-08, TEST_TOL5);
  s += test_sf_check_result(message_buff,  G,  1.333127992006686320e+07, TEST_SQRT_TOL0);
  s += test_sf_check_result(message_buff, Gp, -4.715914530842402330e+07, TEST_SQRT_TOL0);
  printf("%s", message_buff);
  gsl_test(s, "  gsl_sf_coulomb_wave_FG_e(8.0, 1.05, lam_F=0, lam_G=0)");
  status += s;

  lam_F = 0.1;
  k_G = 0;
  eta = 8.0;
  x = 1.05;
  gsl_sf_coulomb_wave_FG_e(eta, x, lam_F, k_G, &F, &Fp, &G, &Gp, &Fe, &Ge);
  s = 0;
  message_buff[0] = 0;
  s += test_sf_check_result(message_buff,  F,  9.611416736061987761e-09, TEST_TOL5);
  s += test_sf_check_result(message_buff, Fp,  3.909628126126824140e-08, TEST_TOL5);
  s += test_sf_check_result(message_buff,  G,  1.365928464219262581e+07, 4.0*TEST_SQRT_TOL0);
  s += test_sf_check_result(message_buff, Gp, -4.848117385783386850e+07, 4.0*TEST_SQRT_TOL0);
  printf("%s", message_buff);
  gsl_test(s, "  gsl_sf_coulomb_wave_FG_e(8.0, 1.05, lam_F=0.1, lam_G=0.1)");
  status += s;

  lam_F = 0.0;
  k_G = 0;
  eta = 50.0;
  x = 0.1;
  gsl_sf_coulomb_wave_FG_e(eta, x, lam_F, k_G, &F, &Fp, &G, &Gp, &Fe, &Ge);
  s = 0;
  message_buff[0] = 0;
  s += test_sf_check_result(message_buff,  F,  2.807788027954216071e-67, TEST_TOL5);
  s += test_sf_check_result(message_buff, Fp,  9.677600748751576606e-66, TEST_TOL5);
  s += test_sf_check_result(message_buff,  G,  5.579810686998358766e+64, TEST_SQRT_TOL0);
  s += test_sf_check_result(message_buff, Gp, -1.638329512756321424e+66, TEST_SQRT_TOL0);
  printf("%s", message_buff);
  gsl_test(s, "  gsl_sf_coulomb_wave_FG_e(50.0, 0.1, lam_F=0, lam_G=0)");
  status += s;

  lam_F = 0.0;
  k_G = 0;
  eta = 10.0;
  x = 5.0;
  gsl_sf_coulomb_wave_FG_e(eta, x, lam_F, k_G, &F, &Fp, &G, &Gp, &Fe, &Ge);
  s = 0;
  message_buff[0] = 0;
  s += test_sf_check_result(message_buff,  F,  1.7207454091787930614e-06, 10.0*WKB_TOL);
  s += test_sf_check_result(message_buff, Fp,  3.0975994706405458046e-06, 10.0*WKB_TOL);
  s += test_sf_check_result(message_buff,  G,  167637.56609459967623, 10.0*WKB_TOL);
  s += test_sf_check_result(message_buff, Gp, -279370.76655361803075, 10.0*WKB_TOL);
  printf("%s", message_buff);
  gsl_test(s, "  gsl_sf_coulomb_wave_FG_e(10.0, 5.0, lam_F=0, lam_G=0)");
  status += s;

  lam_F = 0.0;
  k_G = 0;
  eta = 25.0;
  x = 10.0;
  gsl_sf_coulomb_wave_FG_e(eta, x, lam_F, k_G, &F, &Fp, &G, &Gp, &Fe, &Ge);
  s = 0;
  message_buff[0] = 0;
  s += test_sf_check_result(message_buff,  F,  1.5451274501076114315e-16, 5.0*WKB_TOL);
  s += test_sf_check_result(message_buff, Fp,  3.1390869393378630928e-16, 5.0*WKB_TOL);
  s += test_sf_check_result(message_buff,  G,  1.6177129008336318136e+15, 5.0*WKB_TOL);
  s += test_sf_check_result(message_buff, Gp, -3.1854062013149740860e+15, 5.0*WKB_TOL);
  printf("%s", message_buff);
  gsl_test(s, "  gsl_sf_coulomb_wave_FG_e(25.0, 10.0, lam_F=0, lam_G=0)");
  status += s;

  lam_F = 0.0;
  k_G = 0;
  eta = 1.0;
  x = 9.2;
  gsl_sf_coulomb_wave_FG_e(eta, x, lam_F, k_G, &F, &Fp, &G, &Gp, &Fe, &Ge);
  s = 0;
  message_buff[0] = 0;
  s += test_sf_check_result(message_buff,  F, -0.25632012319757955655, TEST_TOL5);
  s += test_sf_check_result(message_buff, Fp,  0.91518792286724220370, TEST_TOL5);
  s += test_sf_check_result(message_buff,  G,  1.03120585918973466110, TEST_SQRT_TOL0);
  s += test_sf_check_result(message_buff, Gp,  0.21946326717491250193, TEST_SQRT_TOL0);
  printf("%s", message_buff);
  gsl_test(s, "  gsl_sf_coulomb_wave_FG_e(1.0, 9.2, lam_F=0, lam_G=0)");
  status += s;

  lam_F = 0.0;
  eta = 10.0;
  x = 10.0;
  gsl_sf_coulomb_wave_FG_e(eta, x, lam_F, k_G, &F, &Fp, &G, &Gp, &Fe, &Ge);
  s = 0;
  message_buff[0] = 0;
  s += test_sf_check_result(message_buff,  F,  0.0016262711250135878249, WKB_TOL);
  s += test_sf_check_result(message_buff, Fp,  0.0017060476320792806014, WKB_TOL);
  s += test_sf_check_result(message_buff,  G,  307.87321661090837987, WKB_TOL);
  s += test_sf_check_result(message_buff, Gp, -291.92772380826822871, WKB_TOL);
  printf("%s", message_buff);
  gsl_test(s, "  gsl_sf_coulomb_wave_FG_e(10.0, 10.0, lam_F=0, lam_G=0)");
  status += s;

  lam_F = 0.0;
  eta = 100.0;
  x = 1.0;
  gsl_sf_coulomb_wave_FG_e(eta, x, lam_F, k_G, &F, &Fp, &G, &Gp, &Fe, &Ge);
  s = 0;
  message_buff[0] = 0;
  s += test_sf_check_result(message_buff,  F,  8.999367996930662705e-126, 10.0*WKB_TOL);
  s += test_sf_check_result(message_buff, Fp,  1.292746745757069321e-124, 10.0*WKB_TOL);
  s += test_sf_check_result(message_buff,  G,  3.936654148133683610e+123, 10.0*WKB_TOL);
  s += test_sf_check_result(message_buff, Gp, -5.456942268061526371e+124, 10.0*WKB_TOL);
  printf("%s", message_buff);
  gsl_test(s, "  gsl_sf_coulomb_wave_FG_e(100.0, 1.0, lam_F=0, lam_G=0)");
  status += s;

  /* compute F_1(eta=0,x=3.25), F'_1 and G_1(eta=0,x=3.25), G'_1 */

  lam_F = 1.0;
  eta = 0.0;
  x = 3.25;
  k_G = 0;
  gsl_sf_coulomb_wave_FG_e(eta, x, lam_F, k_G, &F, &Fp, &G, &Gp, &Fe, &Ge);
  s = 0;
  message_buff[0] = 0;
  s += test_sf_check_result(message_buff,  F,  sin(x)/x - cos(x), TEST_TOL3);
  s += test_sf_check_result(message_buff, Fp,  -sin(x)/(x*x) + cos(x)/x +sin(x), TEST_TOL3);
  s += test_sf_check_result(message_buff,  G,  cos(x)/x + sin(x), TEST_TOL3);
  s += test_sf_check_result(message_buff, Gp,  -cos(x)/(x*x) - sin(x)/x + cos(x), TEST_TOL3);
  printf("%s", message_buff);
  gsl_test(s, "  gsl_sf_coulomb_wave_FG_e(3.25, 0.0, lam_F=1, lam_G=1)");
  status += s;

  /* compute F_1(eta=0,x=3.25), F'_1 and G_0(eta=0,x=3.25), G'_0 */

  lam_F = 1.0;
  eta = 0.0;
  x = 3.25;
  k_G = 1;
  gsl_sf_coulomb_wave_FG_e(eta, x, lam_F, k_G, &F, &Fp, &G, &Gp, &Fe, &Ge);
  s = 0;
  message_buff[0] = 0;
  s += test_sf_check_result(message_buff,  F,  sin(x)/x - cos(x), TEST_TOL3);
  s += test_sf_check_result(message_buff, Fp,  -sin(x)/(x*x) + cos(x)/x +sin(x), TEST_TOL3);
  s += test_sf_check_result(message_buff,  G,  cos(x), TEST_TOL3);
  s += test_sf_check_result(message_buff, Gp,  -sin(x), TEST_TOL3);
  printf("%s", message_buff);
  gsl_test(s, "  gsl_sf_coulomb_wave_FG_e(3.25, 0.0, lam_F=1, lam_G=0)");
  status += s;

  return status;
}


text:
//A one-stop wrapper for Quadrature calculations
//https://scicomp.stackexchange.com/questions/20786/c-library-for-numerical-intergration-quadrature
#ifndef __FIXEDPOINT_QUADRATURE_H__
#define __FIXEDPOINT_QUADRATURE_H__

#include <iostream>
#include <cmath>

#include <functional>
#include <memory>
#include <utility>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_integration.h>

template < typename F >
class gsl_fixedpoint_quad
{
    protected:
        F f;
        int limit;
        std::unique_ptr < gsl_integration_glfixed_table,
                        std::function < void(gsl_integration_glfixed_table*) >
                        > workspace;

        static double gsl_wrapper(double x, void * p)
        {
            gsl_fixedpoint_quad * t = reinterpret_cast<gsl_fixedpoint_quad*>(p);
            return t->f(x);
        }

    public:
        gsl_fixedpoint_quad(F f, int limit)
            : f(f)
            , limit(limit)
            , workspace(gsl_integration_glfixed_table_alloc(limit), gsl_integration_glfixed_table_free)
        {}
        virtual ~gsl_fixedpoint_quad(){};

        double integrate(double min, double max)
        {
            gsl_function gsl_f;
            gsl_f.function = &gsl_wrapper;
            gsl_f.params = this;

            return gsl_integration_glfixed ( &gsl_f, min, max, workspace.get());
        }
};

template < typename F >
double GSLFixedPointQuadrature(F func,
            std::pair<double,double> const& range,
            int limit = 5)
{
    return gsl_fixedpoint_quad<F>(func, limit).integrate(range.first, range.second);
}

#endif //__FIXEDPOINT_QUADRATURE_H__

ext: h

max_stars_repo_name: muhammadhasyim/pyglasstools

max_issues_repo_name: muhammadhasyim/pyglasstools

max_forks_repo_name: muhammadhasyim/pyglasstools


text:
/*
 * Copyright (c) 2018-2021 Aleksas Mazeliauskas, Stefan Floerchinger,
 *                    Eduardo Grossi, and Derek Teaney
 * All rights reserved.
 *
 * FastReso is distributed under MIT license;
 * see the LICENSE file that should be present in the root
 * of the source distribution, or alternately available at:
 * https://github.com/amazeliauskas/FastReso/
 */
#ifndef FASTRESO_TFastReso_h
#define FASTRESO_TFastReso_h

#include <fstream>
#include <memory>
#include <map>
#include <vector>
#include <string>
#include <gsl/gsl_integration.h>

// The class storing the properties of particles and their irreducible distrubtion
// function components f_i.  Printing out routines are also defined there.
#include "TParticle.h"

// List of parameters controlling the momentum discretization of irreducible components
// See also TParticle class definition
#include "grid_params.h"

//! TFastReso class is the top wrapper class of the fast freeze-out implementation.
class TFastReso {
  public:
    int fParticleNumber;
    //! reads the inputfile of resonances and creates TParticle entry in fParticleData vector.
    //! This routine should be edited by the user to match their input file format.
    //! Currently accepted format is (sstr is a string stream from the input file) is the one used
    //! in particle.data file in THERMINATOR 2 	[arXiv:1102.0273]
    //!  Name   Mass   Gamma   Sppin  Isospin  I3    Nq    Ns    Naq    Nas    Nc    Nac    PDGCode;
    //!  Lines starting with '#' are ignored.
    virtual void read_particles_data(std::string inputname, std::string comps, bool verbose=false) = 0;

    //! The function call to initialized the particle distribution functions in fParticleData. Input parameters
    //! are the freeze-out temperature, chemical potential and (optionally), speed of sound (only needed for bulk perturbations).
    virtual void do_thermal(double Tfo, double MuB=0, double MuI3=0, double MuS=0, double MuC=0, double Cs2=0.14) = 0;

    //! Read in the input file of resonance decays and perfom the decays on the fly. It is important that decays
    //! are mass ordered to included all feeddown contributions to lower mass resonances.
    //! The table format is the same as decays.data in THERMINATOR 2 	[arXiv:1102.0273]
    //! Father Child1 Child2 (optional Child3) BranchingRatio ClebschGordanCoeff
    //! Note that ClebschGordanCoeff term is 0 or 1 identifying if the branching ration needs
    //! to be reweighted by actual Clebsch Gordon coefficients.
    virtual void do_decays(std::string inputname, bool verbose=false) = 0;
  public:
    //TFastReso(std::string comps);
    //! routine to get the pointer to TParticle class for a particle particle name"
    int getNumberOfParticles(){return fParticleNumber;};
    //virtual TParticle * getParticle(std::string name);
    //virtual TParticle * getParticle(int i);
};

// Speed of sound squared as a function of temperature T=x in GeV
inline double cs2FluiduM(double x) {
  return (x*x*(0.0010461910330715387 + x*(-0.016997351850612453 + x*(0.12595528994069613 + (-0.510477729039857 + x)*x)))*
     (3.831255349484765e-8 + x*(0.00002140058802480486 + x*(-0.00008693192041426982 +
             x*(-0.018444802527704172 + x*(0.5368632648594895 + x*(-7.4522565817870365 + x*(62.30776783074574 + x*(-336.8744183413884 + x*(1191.026405027395 + x*(-2590.1225767072065 + 2712.1934046206975*x)))))))))
       ))/(7.313548422127797e-15 + x*(8.356306148115294e-12 + x*(2.3387986355395423e-9 +
           x*(-6.605395132516945e-9 + x*(-1.790148526211301e-6 + x*(8.958041058260661e-6 +
                    x*(0.0014381954623440515 + x*(-0.045458133957877116 + x*(0.726263758547687 +
                             x*(-7.564793448616153 + x*(56.12131032495715 + x*(-308.4534327185232 +
                                      x*(1271.3766801470479 + x*(-3884.3701789306033 + x*(8436.311867208919 + x*(-11805.768858200583 + 8136.580213862093*x))))))))))))))));
}

#endif

max_stars_repo_name: amazeliauskas/FastReso

max_issues_repo_name: amazeliauskas/FastReso

max_forks_repo_name: amazeliauskas/FastReso


text:
/**
 *
 * @file core_zpltmg_chebvand.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @precisions normal z -> c d s
 *
 **/
#include <lapacke.h>
#include <math.h>
#include "common.h"

/***************************************************************************//**
 *
 * @ingroup CORE_PLASMA_Complex64_t
 *
 *  CORE_zpltmg_chebvand is a kernel used in Vandermonde-like matrix generation
 *
 *  See http://www.mathworks.fr/fr/help/matlab/ref/gallery.html#f84-999859
 *
 *  Vandermonde-like matrix for the Chebyshev polynomials
 *
 *  Produces the (primal) Chebyshev Vandermonde matrix based on the vector of
 *  points p, which define where the Chebyshev polynomial is calculated.
 *
 *  If seed != 0, C(i,j) = Ti – 1(p(j)) where Ti – 1 is the Chebyshev
 *  polynomial of degree i – 1, and p is a vector of N equally spaced points on
 *  the interval [0,1].
 *
 *******************************************************************************
 *
 * @param[in] M
 *         The number of rows of the tile A to initialize. M >= 2.
 *
 * @param[in] N
 *         The number of columns of the tile A to initialize. N >= 0.
 *
 * @param[out] A
 *         On entry, the M-by-N tile to be initialized.
 *         On exit, each element of A is defined by:
 *         A(i,j) = Ti – 1(p(j)) where Ti – 1 is the Chebyshev polynomial of
 *         degree i – 1
 *
 * @param[in] LDA
 *         The leading dimension of the tile A. LDA >= max(1,M).
 *
 * @param[in] gN
 *         The global number of columns of the full matrix, A is belonging to. gN >= (n0+gN).
 *
 * @param[in] m0
 *         The index of the first row of tile A in the full matrix. m0 >= 0.
 *
 * @param[in] n0
 *         The index of the first column of tile A in the full matrix. n0 >= 0.
 *
 * @param[in] gN
 *         The global number of columns of the full matrix, A is belonging to. gN >= (n0+gN).
 *
 * @param[in,out] W
 *          Workspace of size 2-by-N, that contains the N triplets:
 *               ( A( m0-2, j), A(m0-1, j) )
 *          On entry, if m == 0, W is uinitialized, otherwise contains the data
 *          described above.
 *          On exit, contains the triplets ( A(m0+M-2, j), A(m0+M-1, j) )
 *
 *******************************************************************************
 *
 * @return
 *         \retval PLASMA_SUCCESS successful exit
 *         \retval <0 if INFO = -k, the k-th argument had an illegal value
 *
 ******************************************************************************/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_zpltmg_chebvand = PCORE_zpltmg_chebvand
#define CORE_zpltmg_chebvand PCORE_zpltmg_chebvand
#endif
int CORE_zpltmg_chebvand( int M, int N, PLASMA_Complex64_t *A, int LDA,
                          int gN, int m0, int n0,
                          PLASMA_Complex64_t *W )
{
    PLASMA_Complex64_t step;
    int i, j, jj;

    /* Check input arguments */
    if (M < 0) {
        coreblas_error(1, "Illegal value of M");
        return -1;
    }
    if (N < 0) {
        coreblas_error(2, "Illegal value of N");
        return -2;
    }
    if ((LDA < max(1,M)) && (M > 0)) {
        coreblas_error(4, "Illegal value of LDA");
        return -4;
    }
    if (m0 < 0) {
        coreblas_error(6, "Illegal value of m0");
        return -6;
    }
    if (n0 < 0) {
        coreblas_error(7, "Illegal value of n0");
        return -7;
    }
    if (gN < n0+N) {
        coreblas_error(5, "Illegal value of gN");
        return -5;
    }

    step = (PLASMA_Complex64_t)1. / (gN - 1.);

    /* Initialize W if required */
    if (m0 == 0) {
        for (j=0, jj=n0; j<N; j++, jj++) {
            W[2*j  ] = 1.;
            W[2*j+1] = jj * step;
        }

        if ( M == 1 ) {
            LAPACKE_zlacpy_work( LAPACK_COL_MAJOR, 'A',
                                 1, N, W, 2, A, LDA );
            return PLASMA_SUCCESS;
        }

        LAPACKE_zlacpy_work( LAPACK_COL_MAJOR, 'A',
                             2, N, W, 2, A, LDA );

        M -= 2;
        A += 2;
    }

    /* Case NB=1, W contains row 0 and 1 and M should be 1 */
    if (m0 == 1) {
        if (M != 1) {
            coreblas_error(1, "Illegal value of M for m0 = 1");
            return -1;
        }

        LAPACKE_zlacpy_work( LAPACK_COL_MAJOR, 'A',
                             1, N, W+1, 2, A, LDA );
        return PLASMA_SUCCESS;
    }

    for (j=0, jj=n0; j<N; j++, jj++) {
        /* First line */
        if (M > 0) {
            A[LDA*j] = 2. * jj * step * W[j*2 + 1] - W[j*2  ];
        }

        /* Second line */
        if (M > 1) {
            A[LDA*j+1] = 2. * jj * step * A[LDA*j  ] - W[j*2+1];
        }

        for (i=2; i<M; i++) {
            A[LDA*j+i] = 2. * jj * step * A[LDA*j+i-1] - A[LDA*j+i-2];
        }
    }

    if ( M == 1 ) {
        cblas_zcopy( N, W+1, 2, W, 2 );
        cblas_zcopy( N, A+M-1, LDA, W+1, 2 );
    }
    else {
        LAPACKE_zlacpy_work( LAPACK_COL_MAJOR, 'A',
                             2, N, A + M - 2, LDA, W, 2 );
    }

    return PLASMA_SUCCESS;
}


text:
#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_math.h>

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

/*! \file b_from_rot_a.c
 *  \brief calculates ror(b) to aloow an initial vector potential
 *
 *  This file is only called once after reading the initial condition.
 *  It is used to calculate b_ini from a initial vector potential, if needed.
 */

#ifdef BFROMROTA

void rot_a(void)
{
  long long ntot, ntotleft;
  int ndone;
  int *noffset, *nbuffer, *nsend, *nsend_local, *numlist, *ndonelist;
  int i, j, n;
  int npleft;
  int maxfill, source;
  int level, ngrp, sendTask, recvTask;
  int place, nexport;
  double dmax1, dmax2, fac;
  double sumt, sumcomm;

  MPI_Status status;

  noffset = (int *) mymalloc(sizeof(int) * NTask);	/* offsets of bunches in common list */
  nbuffer = (int *) mymalloc(sizeof(int) * NTask);
  nsend_local = (int *) mymalloc(sizeof(int) * NTask);
  nsend = (int *) mymalloc(sizeof(int) * NTask * NTask);
  ndonelist = (int *) mymalloc(sizeof(int) * NTask);

  NumSphUpdate = N_gas;

  numlist = (int *) mymalloc(NTask * sizeof(int) * NTask);
  MPI_Allgather(&NumSphUpdate, 1, MPI_INT, numlist, 1, MPI_INT, MPI_COMM_WORLD);
  for(i = 0, ntot = 0; i < NTask; i++)
    ntot += numlist[i];
  myfree(numlist);

  i = 0;			/* beginn with this index */
  ntotleft = ntot;		/* particles left for all tasks together */

  while(ntotleft > 0)
    {
      for(j = 0; j < NTask; j++)
	nsend_local[j] = 0;

      /* do local particles and prepare export list */
      for(nexport = 0, ndone = 0; i < N_gas && nexport < All.BunchSizeDensity - NTask; i++)
	{
	  ndone++;

	  for(j = 0; j < NTask; j++)
	    Exportflag[j] = 0;

	  rot_a_evaluate(i, 0);

	  for(j = 0; j < NTask; j++)
	    {
	      if(Exportflag[j])
		{
		  DensDataIn[nexport].Pos[0] = P[i].Pos[0];
		  DensDataIn[nexport].Pos[1] = P[i].Pos[1];
		  DensDataIn[nexport].Pos[2] = P[i].Pos[2];
		  /* using velocity structure for BPred ... */
		  DensDataIn[nexport].Vel[0] = SphP[i].BPred[0];
		  DensDataIn[nexport].Vel[1] = SphP[i].BPred[1];
		  DensDataIn[nexport].Vel[2] = SphP[i].BPred[2];

		  DensDataIn[nexport].Hsml = PPP[i].Hsml;
		  DensDataIn[nexport].Index = i;
		  DensDataIn[nexport].Task = j;
		  nexport++;
		  nsend_local[j]++;
		}
	    }
	}

      qsort(DensDataIn, nexport, sizeof(struct densdata_in), dens_compare_key);

      for(j = 1, noffset[0] = 0; j < NTask; j++)
	noffset[j] = noffset[j - 1] + nsend_local[j - 1];

      MPI_Allgather(nsend_local, NTask, MPI_INT, nsend, NTask, MPI_INT, MPI_COMM_WORLD);

      /* now do the particles that need to be exported */

      for(level = 1; level < (1 << PTask); level++)
	{
	  for(j = 0; j < NTask; j++)
	    nbuffer[j] = 0;
	  for(ngrp = level; ngrp < (1 << PTask); ngrp++)
	    {
	      maxfill = 0;
	      for(j = 0; j < NTask; j++)
		{
		  if((j ^ ngrp) < NTask)
		    if(maxfill < nbuffer[j] + nsend[(j ^ ngrp) * NTask + j])
		      maxfill = nbuffer[j] + nsend[(j ^ ngrp) * NTask + j];
		}
	      if(maxfill >= All.BunchSizeDensity)
		break;

	      sendTask = ThisTask;
	      recvTask = ThisTask ^ ngrp;

	      if(recvTask < NTask)
		{
		  if(nsend[ThisTask * NTask + recvTask] > 0 || nsend[recvTask * NTask + ThisTask] > 0)
		    {
		      /* get the particles */
		      MPI_Sendrecv(&DensDataIn[noffset[recvTask]],
				   nsend_local[recvTask] * sizeof(struct densdata_in), MPI_BYTE,
				   recvTask, TAG_DENS_A,
				   &DensDataGet[nbuffer[ThisTask]],
				   nsend[recvTask * NTask + ThisTask] * sizeof(struct densdata_in),
				   MPI_BYTE, recvTask, TAG_DENS_A, MPI_COMM_WORLD, &status);
		    }
		}

	      for(j = 0; j < NTask; j++)
		if((j ^ ngrp) < NTask)
		  nbuffer[j] += nsend[(j ^ ngrp) * NTask + j];
	    }

	  for(j = 0; j < nbuffer[ThisTask]; j++)
	    {
	      rot_a_evaluate(j, 1);
	    }

	  MPI_Barrier(MPI_COMM_WORLD);

	  for(j = 0; j < NTask; j++)
	    nbuffer[j] = 0;
	  for(ngrp = level; ngrp < (1 << PTask); ngrp++)
	    {
	      maxfill = 0;
	      for(j = 0; j < NTask; j++)
		{
		  if((j ^ ngrp) < NTask)
		    if(maxfill < nbuffer[j] + nsend[(j ^ ngrp) * NTask + j])
		      maxfill = nbuffer[j] + nsend[(j ^ ngrp) * NTask + j];
		}
	      if(maxfill >= All.BunchSizeDensity)
		break;

	      sendTask = ThisTask;
	      recvTask = ThisTask ^ ngrp;

	      if(recvTask < NTask)
		{
		  if(nsend[ThisTask * NTask + recvTask] > 0 || nsend[recvTask * NTask + ThisTask] > 0)
		    {
		      /* send the results */
		      MPI_Sendrecv(&DensDataResult[nbuffer[ThisTask]],
				   nsend[recvTask * NTask + ThisTask] * sizeof(struct densdata_out),
				   MPI_BYTE, recvTask, TAG_DENS_B,
				   &DensDataPartialResult[noffset[recvTask]],
				   nsend_local[recvTask] * sizeof(struct densdata_out),
				   MPI_BYTE, recvTask, TAG_DENS_B, MPI_COMM_WORLD, &status);

		      /* add the result to the particles */
		      for(j = 0; j < nsend_local[recvTask]; j++)
			{
			  source = j + noffset[recvTask];
			  place = DensDataIn[source].Index;

			  SphP[place].BSmooth[0] += DensDataPartialResult[source].BSmooth[0];
			  SphP[place].BSmooth[1] += DensDataPartialResult[source].BSmooth[1];
			  SphP[place].BSmooth[2] += DensDataPartialResult[source].BSmooth[2];
			  SphP[place].DensityNorm += DensDataPartialResult[source].DensityNorm;
			}
		    }
		}

	      for(j = 0; j < NTask; j++)
		if((j ^ ngrp) < NTask)
		  nbuffer[j] += nsend[(j ^ ngrp) * NTask + j];
	    }

	  level = ngrp - 1;
	}

      MPI_Allgather(&ndone, 1, MPI_INT, ndonelist, 1, MPI_INT, MPI_COMM_WORLD);
      for(j = 0; j < NTask; j++)
	ntotleft -= ndonelist[j];
    }

  for(i = 0; i < N_gas; i++)
    for(j = 0; j < 3; j++)
#ifndef IGNORE_PERIODIC_IN_ROTA
      SphP[i].B[j] = SphP[i].BPred[j] = SphP[i].BSmooth[j] / SphP[i].d.Density;
#else
      SphP[i].B[j] = SphP[i].BPred[j] = SphP[i].BSmooth[j] / SphP[i].DensityNorm;
#endif
  myfree(ndonelist);
  myfree(nsend);
  myfree(nsend_local);
  myfree(nbuffer);
  myfree(noffset);

}

/*! This function represents the core of the calculus of rot of vector potential
 *  target particle may either be local, or reside in the communication
 *  buffer.
 */
void rot_a_evaluate(int target, int mode)
{
  int j, n;
  int startnode, numngb, numngb_inbox;
  double h, h2, fac, hinv, hinv3, hinv4;

  double wk, dwk;
  double dx, dy, dz, r, r2, u, mass_j, rho;
  double dbx, dby, dbz;

  double rotb[3];

  MyFloat *pos, *b;

  rho = rotb[0] = rotb[1] = rotb[2] = 0;

  if(mode == 0)
    {
      pos = P[target].Pos;
      b = SphP[target].BPred;
      h = PPP[target].Hsml;
    }
  else
    {
      pos = DensDataGet[target].Pos;
      b = DensDataGet[target].Vel;	/* Using Vel structure for BPred !!!! */
      h = DensDataGet[target].Hsml;
    }

  h2 = h * h;
  hinv = 1.0 / h;
#ifndef  TWODIMS
  hinv3 = hinv * hinv * hinv;
#else
  hinv3 = hinv * hinv / boxSize_Z;
#endif
  hinv4 = hinv3 * hinv;

  startnode = All.MaxPart;
  do
    {
      numngb_inbox = ngb_treefind_variable(&pos[0], 2 * h, &startnode);

      for(n = 0; n < numngb_inbox; n++)
	{
	  j = Ngblist[n];

	  dx = pos[0] - P[j].Pos[0];
	  dy = pos[1] - P[j].Pos[1];
	  dz = pos[2] - P[j].Pos[2];

	  dbx = b[0] - SphP[j].BPred[0];
	  dby = b[1] - SphP[j].BPred[1];
	  dbz = b[2] - SphP[j].BPred[2];

#ifndef IGNORE_PERIODIC_IN_ROTA
#ifdef PERIODIC			/*  now find the closest image in the given box size  */
	  if(dx > boxHalf_X)
	    dx -= boxSize_X;
	  if(dx < -boxHalf_X)
	    dx += boxSize_X;
	  if(dy > boxHalf_Y)
	    {
#ifdef BRIOWU
	      dbz /= dy;
#endif
	      dy -= boxSize_Y;
#ifdef BRIOWU
	      dbz *= dy;
#endif
	    }
	  if(dy < -boxHalf_Y)
	    {
#ifdef BRIOWU
	      dbz /= dy;
#endif
	      dy += boxSize_Y;
#ifdef BRIOWU
	      dbz *= dy;
#endif
	    }
	  if(dz > boxHalf_Z)
	    {
#ifdef BRIOWU
	      dbx /= dz;
#endif
	      dz -= boxSize_Z;
#ifdef BRIOWU
	      dbx *= dz;
#endif
	    }
	  if(dz < -boxHalf_Z)
	    {
#ifdef BRIOWU
	      dbx /= dz;
#endif
	      dz += boxSize_Z;
#ifdef BRIOWU
	      dbx *= dz;
#endif
	    }
#endif
#endif
	  r2 = dx * dx + dy * dy + dz * dz;

	  if(r2 < h2)
	    {
	      r = sqrt(r2);
	      u = r * hinv;
	      if(u < 0.5)
		{
		  wk = hinv3 * (KERNEL_COEFF_1 + KERNEL_COEFF_2 * (u - 1) * u * u);
		  dwk = hinv4 * u * (KERNEL_COEFF_3 * u - KERNEL_COEFF_4);
		}
	      else
		{
		  wk = hinv3 * KERNEL_COEFF_5 * (1.0 - u) * (1.0 - u) * (1.0 - u);
		  dwk = hinv4 * KERNEL_COEFF_6 * (1.0 - u) * (1.0 - u);
		}

	      mass_j = P[j].Mass;

	      rho += mass_j * wk;

	      if(r > 0)
		{
		  fac = mass_j * dwk / r;
#ifndef BRIOWU
		  rotb[0] += FLT(fac * (dz * dby - dy * dbz));
		  rotb[1] += FLT(fac * (dx * dbz - dz * dbx));
		  rotb[2] += FLT(fac * (dy * dbx - dx * dby));
#else
		  rotb[0] += FLT(fac * (-dy * dbz));
		  rotb[1] += FLT(fac * (-dz * dbx));
		  rotb[2] += FLT(fac * (0.0));
#endif
		}
	    }
	}

    }
  while(startnode >= 0);

  if(mode == 0)
    {
      SphP[target].BSmooth[0] = rotb[0];	/* giving back rot(B) in the Smooth */
      SphP[target].BSmooth[1] = rotb[1];	/* Data structure !!! */
      SphP[target].BSmooth[2] = rotb[2];
      SphP[target].DensityNorm = rho;
    }
  else
    {
      DensDataResult[target].BSmooth[0] = rotb[0];
      DensDataResult[target].BSmooth[1] = rotb[1];
      DensDataResult[target].BSmooth[2] = rotb[2];
      DensDataResult[target].DensityNorm = rho;
    }
}

#endif

max_stars_repo_name: egpbos/egp

max_issues_repo_name: egpbos/egp

max_forks_repo_name: egpbos/egp


text:
/**
* Various binning and weighting routines for aperture pixel tables.
*/

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

#include "aXe_grism.h"
#include "aper_conf.h"
#include "spc_driz.h"
#include "spce_output.h"
#include "spc_optimum.h"
#include "spce_binning.h"

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

#define DEBUG_ME 0x200

#define PIXWEIGHT(x1,y1,x2,y2,pix) ((pix)->weight_function((x1), (y1), (x2),\
  (y2),(pix)))
#define NAIVE_VAL_TO_BININD(x) ((int)floor((x)+1e-6))

/**
 * Function: add_to_spec_table
 * adds some count to an entry in the spectrum table
 *
 * Parameters:
 * @param spec the spectrum table to work on
 * @param bin index of spectrum table entry to add cur_p to
 * @param cur_p the ap_pixel to add
 * @param weight weight of the count
 */
static void
add_to_spec_table (spectrum * const spec, const int bin,
		   const ap_pixel * const cur_p, const int quant_cont,
		   const double weight)
{

  spc_entry *const sp_e = spec->spec + bin;

  // some conditions which should not be violated
  if ((bin < 0) || (bin > spec->spec_len))
    {
      aXe_message (aXe_M_WARN4, __FILE__, __LINE__,
		   "Assignment out of spectrum: %d", bin);
      return;
    }
  if (weight < 0)
    {
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		   "Weight cannot be negative " "but is %f", weight);
    }

  // check whether the spectral element is new
  // and without values up to now
  if (isnan (sp_e->lambda_mean))
    {

      // initialize the spectral element
      sp_e->lambda_mean = cur_p->lambda;
      sp_e->dlambda = cur_p->dlambda;
      sp_e->lambda_max = cur_p->lambda;
      sp_e->lambda_min = cur_p->lambda;
      sp_e->weight = weight;
      sp_e->count = cur_p->count * weight;
      sp_e->error = fabs(cur_p->error) * weight;
      sp_e->dq = cur_p->dq;

      // initialize the contamination
      // this is different for quantitative and
      // geometrical contamination
      if (quant_cont)
	{
	  if ((int)(sp_e->contam==-1)&&((int)cur_p->contam!=-1))
	    {
	      sp_e->contam = cur_p->contam * weight;
	    }
	}
      else
	{
	  if ((int)(sp_e->contam==-1)&&((int)cur_p->contam!=-1))
	    {
	      sp_e->contam = cur_p->contam;
	    }
	}
    }
  else
    {

      // update an existing spectral bin

      // find new maxima and minima
      sp_e->lambda_max = MAX (cur_p->lambda, sp_e->lambda_max);
      sp_e->lambda_min = MIN (cur_p->lambda, sp_e->lambda_min);

      // find new mean lambda via weighted summation
      sp_e->lambda_mean =
	(sp_e->lambda_mean * sp_e->weight +
	 cur_p->lambda * weight) / (weight + sp_e->weight);

      // NEWNEWNEW::
      // find dlambda via weighted summation
      sp_e->dlambda =
	(sp_e->dlambda * sp_e->weight +
	 cur_p->dlambda * weight) / (weight + sp_e->weight);

      // add the weight
      sp_e->weight += weight;

      // add the counts
      sp_e->count += cur_p->count * weight;

      // process the error
      sp_e->error = sqrt (SQR (sp_e->error) + SQR (fabs(cur_p->error) * weight));

      // logically XOR the dq
      sp_e->dq = (sp_e->dq | cur_p->dq);

      // update the contamination,
      // take into account the quantitative
      // contamination
      if (quant_cont)
	{
	  if (((int)sp_e->contam==-1)&&((int)cur_p->contam!=-1))
	    {
	      sp_e->contam = cur_p->contam * weight;
	    }
	  if (((int)sp_e->contam!=-1)&&((int)cur_p->contam!=-1))
	    {
	      sp_e->contam += cur_p->contam * weight;
	    }
	}
      else
	{
	  if (((int)sp_e->contam==-1)&&((int)cur_p->contam!=-1))
	    {
	      sp_e->contam = cur_p->contam;
	    }
	  if (((int)sp_e->contam!=-1)&&((int)cur_p->contam!=-1))
	    {
	      sp_e->contam += cur_p->contam;
	    }
	}
    }
}

/**
 * Function: pixweight_x
 * computes a weight for the bin between coordinates (x1,y1) and (x2,y2)
 * if the pi/4<=angle<=3 pi/4.
 *
 * Parameters:
 * @param x1 - x coordinate for the start point of the bin on the trace
 * @param y1 - y coordinate for the start point of the bin on the trace
 * @param x2 - x coordinate for the end point of the bin on the trace
 * @param y2 - y coordinate for the end point of the bin on the trace
 * @param pix the pixel to compute the weight for in the form of a
 *            w_pixel structure filled out by fill_w_pixel
 *
 * Returns:
 * @return sum - the pixel weight
 */
double
pixweight_x (const double x1, const double y1, const double x2,
	     const double y2, const struct w_pixel_s *const pix)
{
  double a, b;
  double sum = 0;

  a = x1 - pix->cota * (y1 - pix->y0);
  b = x2 - pix->cota * (y2 - pix->y0);
  if ((b >= pix->p0) && (a <= pix->p1))
    {
      sum += pix->slope / 2 * (MIN (b, pix->p1) - MAX (a, pix->p0))
	* (MAX (a, pix->p0) - 2 * pix->p0 + MIN (b, pix->p1))
	* sin (pix->angle);
    }
  if ((b >= pix->p1) && (a <= pix->p2))
    {
      sum += pix->fmax * (MIN (b, pix->p2) - MAX (a, pix->p1))
	* sin (pix->angle);
       }
  if ((b >= pix->p2) && (a <= pix->p3))
    {
      sum += pix->slope / 2 * (MAX (a, pix->p2) - MIN (b, pix->p3))
	* (MAX (a, pix->p2) - 2 * pix->p3 + MIN (b, pix->p3))
	* sin (pix->angle);
    }
  return sum;
}

/**
 * Funtion: pixweight_y
 * Computes a weight for the bin between coordinates (x1,y1) and (x2,y2)
 * if the angle<=pi/4 or 3 pi/4<=angle.
 *
 * @param x1 - x coordinate for the start point of the bin on the trace
 * @param y1 - y coordinate for the start point of the bin on the trace
 * @param x2 - x coordinate for the end point of the bin on the trace
 * @param y2 - y coordinate for the end point of the bin on the trace
 * @param pix - the pixel to compute the weight for in the form of a
 *              w_pixel structure filled out by fill_w_pixel
 *
 * Returns:
 * @return sum - the pixel weight
 */
double
pixweight_y (const double x1, const double y1, const double x2,
	     const double y2, const struct w_pixel_s *const pix)
{
  double a, b;
  double sum = 0;

  a = y1 + pix->tana * (pix->x0 - x1);
  b = y2 + pix->tana * (pix->x0 - x2);
  if (a > b)
    {
      double tmp;
      tmp = a;
      a = b;
      b = tmp;
    }

  if ((b >= pix->p0) && (a <= pix->p1))
    {
      sum += pix->slope / 2 * (MIN (b, pix->p1) - MAX (a, pix->p0))
	* (MAX (a, pix->p0) - 2 * pix->p0 + MIN (b, pix->p1))
	* cos (pix->angle);
    }
  if ((b >= pix->p1) && (a <= pix->p2))
    {
      sum += pix->fmax * (MIN (b, pix->p2) - MAX (a, pix->p1))
	* cos (pix->angle);
    }
  if ((b >= pix->p2) && (a <= pix->p3))
    {
      sum += pix->slope / 2 * (MAX (a, pix->p2) - MIN (b, pix->p3))
	* (MAX (a, pix->p2) - 2 * pix->p3 + MIN (b, pix->p3))
	* cos (pix->angle);
    }
  return sum;
}

/**
 * Function: fill_w_pixel
 * precomputes some properties of a given pixel for purposes
 * of computing the weights it contributes to a given xi bin.
 *
 * Parameters:
 * @param pix a pointer to the w_pix structure to fill
 * @param x0 the x coordinate of the pixel's lower left corner
 * @param y0 the y coordinate of the pixel's lower left corner
 * @param size the size of the pixel
 * @param angle the orientation of the object that has generated the
 *        spectrum
 */
/* since we're only interested in weights here, the size of the square
 *  doesn't matter.  Where I left explicit ones, you could write size
 * and get the true area instead of the fraction of the total area.
 * Dunno why you'd want to do this, though.
 */
static void
fill_w_pixel (w_pixel * const pix, const double x0, const double y0,
	      const double angle)
{
  pix->tana = tan (angle);
  pix->cota = 1 / pix->tana;
  pix->angle = angle;
  pix->x0 = x0;
  pix->y0 = y0;

  if ((angle >= M_PI / 4) && (angle <= 3 * M_PI / 4))
    {
      pix->p0 = MIN (x0, x0 - 1 * pix->cota);
      pix->p1 = MAX (x0, x0 - 1 * pix->cota);
      pix->p2 = MIN (x0 + 1, x0 + 1 - 1 * pix->cota);
      pix->p3 = MAX (x0 + 1, x0 + 1 - 1 * pix->cota);
      pix->fmax = 1 / sin (angle);

      if (fabs (pix->p1 - pix->p0) < 1e-7)
	pix->slope = 0;
      else
	pix->slope = pix->fmax / (pix->p1 - pix->p0);

      pix->weight_function = &pixweight_x;
    }
  else
    {			/* angle between 0 and pi/4 or 3*pi/4 and pi */
      pix->p0 = MIN (y0, y0 - 1 * pix->tana);
      pix->p1 = MAX (y0, y0 - 1 * pix->tana);
      pix->p2 = MIN (y0 + 1, y0 + 1 - 1 * pix->tana);
      pix->p3 = MAX (y0 + 1, y0 + 1 - 1 * pix->tana);
      pix->fmax = 1 / cos (angle);

      if (fabs (pix->p1 - pix->p0) < 1e-7)
	pix->slope = 0;
      else
	pix->slope = pix->fmax / (pix->p1 - pix->p0);

      pix->weight_function = &pixweight_y;
    }
}

/**
 * Function: bin_naive
 * computes a spectrum from a table of aperture pixels generated from
 * spc_extract.  This is adds up the pixel values, distributing them
 * over the trace, taking into account the fracton of the pixel that
 * projects onto the given [xi,xi+1] interval.
 * Return NULL if ap_p is NULL
 *
 * Parameters:
 * @param ap_p the table of aperture pixels
 * @param px_width width of a pixel (=1 if not subsampled)
 * @param ob_width the width of the object that has generated the spectrum
 * @param ob_orientation the orientation of the object that has
 *             generated the spectrum
 * @param flags problems that were accumulated in generating ap_p;
 *   possible flags are defined for the warning member of the spectrum struct
 *
 * Returns:
 * @return spec - the 1D spectrum
 */
spectrum *
bin_naive (const ap_pixel * const ap_p, const double ob_width,
	   const double ob_orient, const int quant_cont)
{
  const ap_pixel *cur_p;
  int upper, lower;

  int bin;
  spectrum *spec, *tspec;
  double phi_trace;
  spc_entry *spec_table;
  double d, frac;

  // immediately return empty PET's
  if (ap_p==NULL)
    return NULL;

  // go through the PET,
  // find the minimum and
  // maximum in trace distance
  cur_p = ap_p;
  upper = NAIVE_VAL_TO_BININD (cur_p->xi);
  lower = NAIVE_VAL_TO_BININD (cur_p->xi);
  while (cur_p->p_x != -1)
    {
      bin = NAIVE_VAL_TO_BININD (cur_p->xi);
      upper = MAX (bin, upper);
      lower = MIN (bin, lower);
      cur_p++;
    }

  // check whether the spectrum
  // will ahve a finite length,
  // exit if not
  if (upper == lower)
    {
      aXe_message (aXe_M_WARN4, __FILE__, __LINE__,
		   "Pixel table empty.\n");
      return NULL;
    }

  // Thresh in some headway so we don't need to worry too much about
  //     accessing invalid elements now and then
  lower -= 10;
  upper += 10;
  spec = allocate_spectrum (upper - lower);

  spec_table = spec->spec;
  cur_p = ap_p;

  while (cur_p->p_x != -1)
    {
      // Compute any fractional pixel that might fall within the
      //desired extraction width
      d = ob_width;
      frac = 1.;

      // continue if the pixel is outside
      // of the extraction region
      if (fabs (cur_p->dist) > d + .5)
	{
	  cur_p++;
	  continue;
	}
      if ((fabs (cur_p->dist) >= d - .5) && (fabs (cur_p->dist) <= d + .5))
	{
	  frac = fabs (d - (fabs (cur_p->dist) - 0.5));
	}

      // store the local trace angle
      phi_trace = cur_p->dxs;

      if (1)
	{
	  double xi;
	  w_pixel pix;
	  double sinp = sin (phi_trace);
	  double cosp = cos (phi_trace);
	  double xc, yc;
	  double w;
	  double sum = 0;

	  fill_w_pixel (&pix, cur_p->x,	cur_p->y, ob_orient);

	  xc = cur_p->xs;
	  yc = cur_p->ys;

	  // at cur_p->xi, there has to be some contribution. We go back
	  // collecting, until w is zero
	  for (xi = cur_p->xi;; xi -= 1)
	    {
	      bin = NAIVE_VAL_TO_BININD (xi);
	      w = PIXWEIGHT (xc + (bin - cur_p->xi) * cosp,
			     yc + (bin - cur_p->xi) * sinp,
			     xc + (bin + 1 - cur_p->xi) * cosp,
			     yc + (bin + 1 - cur_p->xi) * sinp,
			     &pix);
	      if (w < 1e-10)
		break;

	      add_to_spec_table (spec, bin - lower, cur_p, quant_cont,
				 w * frac * cur_p->weight);
	      //	      if (cur_p->weight > 10.0){
	      //		fprintf(stdout,"weight: %f, distance: %f, ewidth: %f\n",
	      //			cur_p->weight, cur_p->dist, d+0.5);
	      //	      }
	      sum += w;
	    }

	  /* Now collect contributions upward of cur_p->xi */
	  for (xi = cur_p->xi + 1;; xi += 1)
	    {
	      bin = NAIVE_VAL_TO_BININD (xi);

	      w = PIXWEIGHT (xc + (bin - cur_p->xi) * cosp,
			     yc + (bin - cur_p->xi) * sinp,
			     xc + (bin + 1 - cur_p->xi) * cosp,
			     yc + (bin + 1 - cur_p->xi) * sinp,
			     &pix);
	      if (w < 1e-10)
		break;

	      add_to_spec_table (spec, bin - lower, cur_p, quant_cont,
				 w * frac * cur_p->weight);
	      sum += w;
	    }

	  if (fabs (sum - 1) > 1e-6)
	    {
	      fprintf(stdout,
		      "Weights added up to only %f for pixel from %4d,%4d\n",
		      sum, cur_p->p_x, cur_p->p_y);
	    }
	}
      cur_p++;
    }

  /* Trimming the INDEF beginning and ending values in spectrum */
  tspec = trim_spectrum (spec);
  free_spectrum (spec);
  spec = NULL;

  // return the spectrum
  return tspec;
}

/**
 * Function: bin_optimal
 * computes a spectrum from a table of aperture pixels generated from
 * spc_extract.  This is adds up the pixel values, distributing them
 * over the trace, taking into account the fracton of the pixel that
 * projects onto the given [xi,xi+1] interval.
 * Return NULL if ap_p is NULL
 *
 * Parameters:
 * @param ap_p the table of aperture pixels
 * @param px_width width of a pixel (=1 if not subsampled)
 * @param ob_width the width of the object that has generated the spectrum
 * @param ob_orientation the orientation of the object that has
 *             generated the spectrum
 * @param flags problems that were accumulated in generating ap_p;
 *   possible flags are defined for the warning member of the spectrum struct
 *
 * Returns:
 * @return spec - the 1D spectrum
 */
spectrum *
bin_optimal (const ap_pixel * const ap_p, const beam curbeam,
	     const int quant_cont, const gsl_matrix *weights,
	     const drzstamp_dim dimension,gsl_matrix *coverage)
{

  const ap_pixel *cur_p;
  ap_pixel *tmp_p;
  spectrum *spec;
  spectrum *tspec;
  spc_entry *spec_table;

  quadrangle quad;
  //  drzstamp_dim dimension;

  double jacob, arr;
  double frac, totweight;

  int jcen, icen;
  int jupp, iupp;
  int jlow, ilow;

  int ii, jj;
  int stpi, stpj;

  // return NULL if
  // empty PET
  if (ap_p==NULL)
    return NULL;

  // allocate memory
  tmp_p = (ap_pixel *) malloc(sizeof(ap_pixel));

  // allocate memory for the spectrum
  spec = allocate_spectrum (weights->size1);
  spec_table = spec->spec;

  // go over each PET pixel
  cur_p = ap_p;
  for (cur_p = ap_p; cur_p->p_x != -1; cur_p++)
    {

      // continue if the pixel is outside
      // of the extraction region
      if (fabs (cur_p->dist) > curbeam.width + .5)
	  continue;

      // determine which fraction
      // of the pixel is inside of the extraction area
      if ((fabs (cur_p->dist) >= curbeam.width - .5) && (fabs (cur_p->dist) <= curbeam.width + .5))
	frac = fabs (curbeam.width - (fabs (cur_p->dist) - 0.5));
      else
	frac = 1.;

      // transfer values to the temporary pixel
      tmp_p->lambda  = cur_p->xi;
      tmp_p->dist    = cur_p->dist;
      tmp_p->dxs     = cur_p->dxs;
      tmp_p->dlambda = 1.0;

      // create the quadrangle for the current pixel
      quad = get_quad_from_pixel(tmp_p, curbeam.orient, dimension);

      // get the jacobian (well, easy here)
      // the term "cos(cur_p->dxs)" must be there
      // to correct the enlargement necessary
      // to cover the whole lambda-crossdispersion area!
      // NOT COMPLETELY understood
      jacob = cos(cur_p->dxs);

      // get the central pixel (icen, jcen) of the current PET-pixel
      icen = (int) floor(cur_p->xi   - dimension.xstart+.5);
      jcen = (int) floor(cur_p->dist - dimension.ystart+.5);

      // get the uper and lower extend of the quadrangle in x
      iupp = (int)floor(quad.xmax - (double)icen + 0.5)+1;
      ilow = (int)floor(quad.xmin - (double)icen + 0.5);

      // get the uper and lower extend of the quadrangle in x
      jupp = (int)floor(quad.ymax - (double)jcen + 0.5)+1;
      jlow = (int)floor(quad.ymin - (double)jcen + 0.5);

      // go over the extend in x
      for (ii=ilow;ii<iupp;ii++)
	{
	  // go over the extend in x
	  for (jj=jlow;jj<jupp;jj++)
	    {

	      // get the coordinates of the current output pixel
	      stpi = icen+ii;
	      stpj = jcen+jj;

	      // check whether the current output pixel is within
	      // the stamp image; continue if not
	      if ( (stpi>=dimension.xsize)||(stpi<0)||(stpj>=dimension.ysize)||(stpj<0) )
		continue;

	      // get the area which falls onto the current output pixel
	      arr = boxer(stpi,stpj,quad.x,quad.y);

	      // compute the pixel weight from
	      // the various inputs
	      totweight =  arr*frac*jacob*gsl_matrix_get(weights, stpi, stpj);
	      //totweight =  arr*frac*jacob;//*gsl_matrix_get(weights, stpi, stpj);

	      //gsl_matrix_set(coverage, stpi, stpj, gsl_matrix_get(coverage, stpi, stpj) + arr*frac*jacob);

	      // add the pixel to the spectrum
	      if (totweight > 0.0)
		add_to_spec_table (spec, stpi, cur_p, quant_cont,totweight);
	      //	      gsl_matrix_set(coverage, stpi, stpj, sqrt (SQR (gsl_matrix_get(coverage, stpi, stpj)) + SQR (fabs(cur_p->error) * totweight)));
	    }
	}
    }

  /* Trimming the INDEF beginning and ending values in spectrum */
  tspec = trim_spectrum (spec);
  free_spectrum (spec);
  spec = NULL;

  free(tmp_p);

  // return the spectrum
  return tspec;
}

/**
  does a straight forward summation/binning of an aperture pixel table
  with appropriate weights (cf. Hornes 1986)

  @param ap_p the table of aperture pixels
  @param ob_orient the orientation of the object that has
    generated the spectrum
  @param flags problems that were accumulated in generating ap_p;
    possible flags are defined for the warning member of the spectrum struct
  @param n_sub subsampling factor
  @return an allocated spectrum
  @see spectrum
*/

spectrum *
bin_weighted (const ap_pixel * const ap_p, const double ob_orient,
	      const trace_func * const trace, const int n_sub,
	      const int flags)
{
     gsl_vector *binned_table;
     gsl_vector *wei_table;
     gsl_vector *wei2_table;

     double xi, d, w, wei, wei2;
     int xii, num_bin, bin;
     const ap_pixel *cur_p = ap_p;
     double upper = cur_p->xi, lower = cur_p->xi;
     spectrum *spec;

     while (cur_p->p_x != -1)
       {
	    upper = MAX (cur_p->xi, upper);
	    lower = MIN (cur_p->xi, lower);
	    cur_p++;
       }

     lower -= 10;
     upper += 10;

     lower = floor (lower);
     upper = floor (upper + 1);

     num_bin = floor ((upper - lower) / n_sub);
     spec = allocate_spectrum (num_bin);

     binned_table = gsl_vector_alloc (num_bin);
     wei_table = gsl_vector_alloc (num_bin);
     wei2_table = gsl_vector_alloc (num_bin);

     gsl_vector_set_all (binned_table, 0);
     gsl_vector_set_all (wei_table, 0);
     gsl_vector_set_all (wei2_table, 0);

     cur_p = ap_p;
     while (cur_p->p_x != -1)
       {
	    xi = (cur_p->xi - lower) / n_sub;
	    xii = floor (xi);

	    d = fabs(cur_p->dist);
	    w = exp (-d * d / (2 * 6.66));
	    w = 1.;
	    add_to_spec_table (spec, xii - 1,cur_p, 0, w * (1 - (xi - xii)));
	    add_to_spec_table (spec, xii   , cur_p, 0, w * (xi - xii));

	    gsl_vector_set (wei_table, xii - 1,
			    gsl_vector_get (wei_table, xii - 1) + w);
	    gsl_vector_set (wei_table, xii,
			    gsl_vector_get (wei_table, xii) + w);

	    gsl_vector_set (wei2_table, xii - 1,
			    gsl_vector_get (wei2_table, xii - 1) + w * w);
	    gsl_vector_set (wei2_table, xii,
			    gsl_vector_get (wei2_table, xii) + w * w);

	    //fprintf(stderr,"%f,%d,%d\n",xi,xii,xii+1);
	    cur_p++;
       }
     cur_p = ap_p;
     for (bin = 0; bin < num_bin; bin++)
       {
	    wei = gsl_vector_get (wei_table, bin);
	    wei2 = gsl_vector_get (wei2_table, bin);
	    if (wei2 != 0)
	      {
		   spec->spec[bin].count = spec->spec[bin].count * wei / wei2;
	      }
       }

     spec->warning = 0;
     return spec;
}

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


text:
#pragma once
#include "nkg_point.h"
#include "nkg_rect.h"
#include "utility/types.h"
#include <stb_image.h>
#include <gsl/gsl>
#include <memory>
#include <vector>
#include <stdlib.h>

namespace cws80 {

class im_image {
public:
    im_image();
    ~im_image();

    void load(const char *filename, uint desired_channels);
    void load_from_memory(const u8 *buffer, uint length, uint desired_channels);
    void load_from_memory(gsl::span<const u8> memory, uint desired_channels)
    {
        load_from_memory(memory.data(), memory.size(), desired_channels);
    }

    const u8 *pixel(uint x, uint y) const;
    u8 alpha(uint x, uint y) const;

    bool is_transparent_column(uint x) const;
    bool is_transparent_row(uint y) const;

    im_image cut(const im_recti &r) const;
    std::vector<im_image> hsplit() const;
    std::vector<im_image> vsplit() const;

    im_recti rect() const;
    im_recti crop_alpha_border() const;

    const u8 *data() const { return data_.get(); }
    uint width() const { return w_; }
    uint height() const { return h_; }
    uint channels() const { return c_; }

    im_image(im_image &&other) = default;
    im_image &operator=(im_image &&other) = default;

private:
    uint w_{}, h_{}, n_{}, c_{};
    std::unique_ptr<u8[], void (*)(void *)> data_{nullptr, &stbi_image_free};

public:
    class exception : public std::runtime_error {
    public:
        using std::runtime_error::runtime_error;
    };
};

}  // namespace cws80


text:
/*
 * author: Achim Gaedke
 * created: January 2001
 * file: pygsl/src/histogrammodule.c
 * $Id: matrixmodule.c,v 1.2 2003/07/25 17:15:16 schnizer Exp $
 */

#include <Python.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <pygsl/error_helpers.h>
/*
 * matrix
 */

/* my typedef */

staticforward PyTypeObject matrix_matrixType;

typedef struct {
    PyObject_HEAD
    gsl_matrix* m;
} matrix_matrixObject;

static PyObject* matrix_matrix_get(PyObject *self,
				   PyObject *args
				   )
{
  gsl_matrix* matrix;
  long int i;
  long int j;
  matrix=(gsl_matrix*)((matrix_matrixObject*)self)->m;
  if (matrix==NULL) {
    PyErr_SetString(PyExc_RuntimeError, "matrix.get got a NULL pointer");
    return NULL;
  }

  /* get index */
  if (0==PyArg_ParseTuple(args,"ll",&i,&j))
      return NULL;
  if (i<0 || i>=matrix->size1) {
    gsl_error ("1st index lies outside valid range of 0 .. size1 - 1",
	       __FILE__,
	       __LINE__,
	       GSL_EDOM );
    return NULL;
  }
  if (j<0 || j>=matrix->size2) {
    gsl_error ("2nd index lies outside valid range of 0 .. size2 - 1",
	       __FILE__,
	       __LINE__,
	       GSL_EDOM );
    return NULL;
  }
  return PyFloat_FromDouble(gsl_matrix_get(matrix,i,j));
}

static PyObject* matrix_matrix_set(PyObject *self,
				   PyObject *args
				   )
{
  gsl_matrix* matrix;
  long int i;
  long int j;
  double x;
  matrix=(gsl_matrix*)((matrix_matrixObject*)self)->m;
  if (matrix==NULL) {
    PyErr_SetString(PyExc_RuntimeError, "matrix.get got a NULL pointer");
    return NULL;
  }

  /* get index */
  if (0==PyArg_ParseTuple(args,"lld",&i,&j,&x))
      return NULL;
  if (i<0 || i>=matrix->size1) {
    gsl_error ("1st index lies outside valid range of 0 .. size1 - 1",
	       __FILE__,
	       __LINE__,
	       GSL_EDOM );
    return NULL;
  }
  if (j<0 || j>=matrix->size2) {
    gsl_error ("2nd index lies outside valid range of 0 .. size2 - 1",
	       __FILE__,
	       __LINE__,
	       GSL_EDOM );
    return NULL;
  }
  gsl_matrix_set(matrix,i,j,x);
  Py_INCREF(Py_None);
  return Py_None;
}

static PyObject* matrix_matrix_set_all(PyObject *self,
				       PyObject* args)
{
  gsl_matrix* matrix;
  double x;
  matrix=(gsl_matrix*)((matrix_matrixObject*)self)->m;
  if (0==PyArg_ParseTuple(args,"d",&x))
      return NULL;
  if (matrix==NULL) {
    PyErr_SetString(PyExc_RuntimeError, "matrix.size got a NULL pointer");
    return NULL;
  }
  gsl_matrix_set_all(matrix,x);
  Py_INCREF(Py_None);
  return Py_None;
}

static PyObject* matrix_matrix_set_identity(PyObject *self)
{
  gsl_matrix* matrix;
  matrix=(gsl_matrix*)((matrix_matrixObject*)self)->m;
  if (matrix==NULL) {
    PyErr_SetString(PyExc_RuntimeError, "matrix.size got a NULL pointer");
    return NULL;
  }
  gsl_matrix_set_identity(matrix);
  Py_INCREF(Py_None);
  return Py_None;
}

static PyObject* matrix_matrix_set_zero(PyObject *self)
{
  gsl_matrix* matrix;
  matrix=(gsl_matrix*)((matrix_matrixObject*)self)->m;
  if (matrix==NULL) {
    PyErr_SetString(PyExc_RuntimeError, "matrix.size got a NULL pointer");
    return NULL;
  }
  gsl_matrix_set_zero(matrix);
  Py_INCREF(Py_None);
  return Py_None;
}

static PyObject* matrix_matrix_transpose(PyObject *self)
{
  gsl_matrix* matrix;
  matrix=(gsl_matrix*)((matrix_matrixObject*)self)->m;
  if (matrix==NULL) {
    PyErr_SetString(PyExc_RuntimeError, "matrix.size got a NULL pointer");
    return NULL;
  }
  gsl_matrix_transpose(matrix);
  Py_INCREF(Py_None);
  return Py_None;
}

static PyObject* matrix_matrix_size(PyObject *self)
{
  gsl_matrix* matrix;

  matrix=(gsl_matrix*)((matrix_matrixObject*)self)->m;
  if (matrix==NULL) {
    PyErr_SetString(PyExc_RuntimeError, "matrix.size got a NULL pointer");
    return NULL;
  }

  return Py_BuildValue("(ll)",matrix->size1,matrix->size2);
}

static PyObject* matrix_matrix_add(PyObject *self,
				   PyObject *arg) {
  gsl_matrix* matrix=(gsl_matrix*)((matrix_matrixObject*)self)->m;
  gsl_matrix* matrix2;
  int result;

  if (matrix==NULL) {
    PyErr_SetString(PyExc_RuntimeError, "matrix.add got a NULL pointer");
    return NULL;
  }
  Py_INCREF(arg);
  /* get histogram from argument */
  if (arg==NULL || !PyObject_TypeCheck(arg,&matrix_matrixType)) {
    PyErr_SetString(PyExc_TypeError, "matrix.add requires pygsl.matrix.matrix");
    Py_DECREF(arg);
    return NULL;
  }
  matrix2=(gsl_matrix*)((matrix_matrixObject*)arg)->m;
  if (matrix2==NULL) {
    PyErr_SetString(PyExc_RuntimeError, "matrix in argument has a NULL pointer");
    return NULL;
  }
  result=gsl_matrix_add(matrix,matrix2);
  Py_DECREF(arg);
  if (result!=GSL_SUCCESS) {
    return NULL;
  }
  Py_INCREF(Py_None);
  return Py_None;
}

static PyMethodDef matrix_matrix_methods[] = {
  /*  {"alloc",(PyCFunction)histogram_histogram_alloc,METH_VARARGS,"allocate necessary space"}, */
  {"get",(PyCFunction)matrix_matrix_get,METH_VARARGS,"gets the (i,j) element"},
  {"size",(PyCFunction)matrix_matrix_size,METH_NOARGS,"returns the size as tuple"},
  {"add",(PyCFunction)matrix_matrix_add,METH_O,"adds the given matrix"},
  {"set",(PyCFunction)matrix_matrix_set,METH_VARARGS,"sets the (i,j) element"},
  {"set_zero",(PyCFunction)matrix_matrix_set_zero,METH_NOARGS,"sets all to 0"},
  {"set_identity",(PyCFunction)matrix_matrix_set_identity,METH_NOARGS,"sets diagonal to 1 and others to 0"},
  {"set_all",(PyCFunction)matrix_matrix_set_identity,METH_O,"sets all entries to x"},
  {"transpose",(PyCFunction)matrix_matrix_transpose,METH_NOARGS,"transposes the matrix in place"},
  {NULL, NULL, 0, NULL}           /* sentinel */
};

static int matrix_matrix_init(PyObject *self,
			      PyObject *args,
			      PyObject *kwds
			      )
{

  PyObject* orig_matrix;
  long int n1;
  long int n2;

  static char *kwlist1[] = {"matrix",NULL};
  static char *kwlist2[] = {"size1","size2",NULL};

  /* initialise histogram pointer */
  ((matrix_matrixObject*)self)->m=NULL;

  /* test two call alternatives: */

  if (PyArg_ParseTupleAndKeywords(args, kwds,"O!",kwlist1,&matrix_matrixType,&orig_matrix)) {
    gsl_matrix* matrix;
    gsl_matrix* matrix2=(gsl_matrix*)((matrix_matrixObject*)orig_matrix)->m;
    if (matrix2==NULL) {
      PyErr_SetString(PyExc_RuntimeError, "matrix in argument has a NULL pointer");
      return -1;
    }
    matrix=gsl_matrix_calloc(matrix2->size1,matrix2->size2);
    if (matrix==NULL)
      return -1;
    gsl_matrix_memcpy(matrix,matrix2);
    ((matrix_matrixObject*)self)->m=matrix;
    return 0;
  }
  else {
    /* reset exception */
    PyErr_Clear();
  }
  if (PyArg_ParseTupleAndKeywords(args, kwds,"ll",kwlist2,&n1,&n2)) {
    gsl_matrix* matrix;
    if (n1<=0 || n2<=0) {
      gsl_error ("matrix length must be positive",
		 __FILE__,
		 __LINE__,
		 GSL_EDOM );
      return -1;
    }
    matrix=gsl_matrix_calloc(n1,n2);
    if (matrix==NULL)
      return -1;
    ((matrix_matrixObject*)self)->m=matrix;
    return 0;

  }
  else {
    /* set own error message */
    PyErr_Clear();
    PyErr_SetString(PyExc_TypeError, "matrix.init requires pygsl.matrix.matirx object or two long int argument");
    return -1;
  }
}

static PyObject*
matrix_matrix_getattr(PyObject* obj, char *name)
{
  return Py_FindMethod(matrix_matrix_methods, obj, name);
}

static void
matrix_matrix_dealloc(PyObject* self)
{
  gsl_matrix* matrix;
  matrix=(gsl_matrix*)((matrix_matrixObject*)self)->m;
  if (matrix!=NULL) gsl_matrix_free(matrix);
  self->ob_type->tp_free(self);
}

static
PyTypeObject matrix_matrixType = {
	PyObject_HEAD_INIT(NULL)
	0,
	"pygsl.matrix.matrix",
	sizeof(matrix_matrixObject),
	0,
	(destructor)matrix_matrix_dealloc,      /* tp_dealloc */
	0,                                      /* tp_print */
	matrix_matrix_getattr,                  /* tp_getattr */
	0,					/* tp_setattr */
	0,			                /* tp_compare */
	0,                  			/* tp_repr */
	0,					/* tp_as_number */
	0,	                		/* tp_as_sequence */
	0,    /* &matrix_marix_as_mapping,*/	/* tp_as_mapping */
	0,				        /* tp_hash */
	0,					/* tp_call */
	0,					/* tp_str */
	0,		                        /* tp_getattro */
	0,					/* tp_setattro */
	0,					/* tp_as_buffer */
	Py_TPFLAGS_DEFAULT,		        /* tp_flags */
        0,				        /* tp_doc */
	0, 		                        /* tp_traverse */
	0,			                /* tp_clear */
	0,              			/* tp_richcompare */
	0,					/* tp_weaklistoffset */
	0,			                /* tp_iter */
	0,					/* tp_iternext */
	0,				        /* tp_methods */
	0,					/* tp_members */
	0,					/* tp_getset */
	0,					/* tp_base */
	0,					/* tp_dict */
	0,					/* tp_descr_get */
	0,					/* tp_descr_set */
	0,					/* tp_dictoffset */
	(initproc)matrix_matrix_init,	        /* tp_init */
	NULL,              			/* tp_alloc */
	NULL,                			/* tp_new */
	NULL         			        /* tp_free */
};

static PyMethodDef matrixMethods[] = {
  {NULL, NULL, 0, NULL}        /* Sentinel */
};

void
initmatrix(void)
{
  PyObject* m;
  m=Py_InitModule("matrix", matrixMethods);
  init_pygsl();
  /* init histogram type */
  matrix_matrixType.ob_type  = &PyType_Type;
  matrix_matrixType.tp_alloc = PyType_GenericAlloc;
  matrix_matrixType.tp_new   = PyType_GenericNew;
  matrix_matrixType.tp_free  = _PyObject_Del;

  /* install histogram type */
  /* important! must increment histogram type reference counter */
  Py_INCREF((PyObject*)&matrix_matrixType);
  PyModule_AddObject(m,"matrix", (PyObject*)&matrix_matrixType);
}

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


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

 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice,
 * this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 * this list of conditions and the following disclaimer in the documentation
 * and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 *
 */

/*
 *
 *
 *       Filename:  red_gen.c
 *
 *    Description:  Redfield-type Generator
 *
 *        Version:  1.0
 *        Created:  25/04/2014 22:28:19
 *       Revision:  none
 *        License:  BSD
 *
 *         Author:  Giuseppe Argentieri (ga), giuseppe.argentieri@ts.infn.it
 *   Organization:  Università degli Studi di Trieste
 *
 *
 */

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_math.h>
#include <math.h>
#include "funcs.h"

/*
 *      FUNCTION
 *         Name:  red_mat
 *  Description:  Creating the Redfield generator matrix in Bloch form, taking as
 *  		  arguments all the integrals, the physical parameters
 *  		  and the pointer to the matrix
 *
 */
int red_mat ( gsl_matrix* red_mx , void* params )
{
	double integrals[12] ;
	if ( integration(params,integrals) != 0 )
		return -1;

	/* Setting the integrals */
	double rcc, icc, rss, iss, rsc, isc, rcs, ics, rc0, ic0, rs0, is0 ;
	rcc = integrals[0] ;
	icc = integrals[1] ;
	rss = integrals[2] ;
	iss = integrals[3] ;
	rsc = integrals[4] ;
	isc = integrals[5] ;
	rcs = integrals[6] ;
	ics = integrals[7] ;
	rc0 = integrals[8] ;
	ic0 = integrals[9] ;
	rs0 = integrals[10] ;
	is0 = integrals[11] ;

	/* Copying the parameters */
	struct f_params* pars = (struct f_params*) params ;
	double o_c, b, O, o_1 ;
	assign_p ( pars, &o_c, &b, &O, &o_1 ) ;
	double D = sqrt(POW_2(o_1)-POW_2(O)) ;

	/* Ensuring that the matrix is zeroed */
	gsl_matrix_set_zero ( red_mx ) ;

	/* Building the matrix */
	gsl_matrix_set ( red_mx, 1, 0, (D/o_1)*(-(O/o_1)*ic0+iss+(O/o_1)*icc) ) ;
	gsl_matrix_set ( red_mx, 1, 1, POW_2(D/o_1)*rc0+(O/o_1)*rss+POW_2(O/o_1)*rcc ) ;
	gsl_matrix_set ( red_mx, 1, 2, o_1/4+(O/o_1)*rsc-POW_2(O/o_1)*rcs ) ;
	gsl_matrix_set ( red_mx, 1, 3, (D/o_1)*(rsc-(O/o_1)*rcs) ) ;
	gsl_matrix_set ( red_mx, 2, 0, (D/o_1)*(is0+isc-(O/o_1)*ics) ) ;
	gsl_matrix_set ( red_mx, 2, 1, -(o_1/4)-(O/o_1)*rsc+rcs ) ;
	gsl_matrix_set ( red_mx, 2, 2, POW_2(D/o_1)*rc0+rcc+(O/o_1)*rss ) ;
	gsl_matrix_set ( red_mx, 2, 3, -(D/o_1)*((O/o_1)*rcc+rss) ) ;
	gsl_matrix_set ( red_mx, 3, 0, (1+POW_2(O/o_1))*ics-2*(O/o_1)*isc ) ;
	gsl_matrix_set ( red_mx, 3, 1, -(D/o_1)*rs0 ) ;
	gsl_matrix_set ( red_mx, 3, 2, -(O*D/(o_1*o_1))*rc0 ) ;
	gsl_matrix_set ( red_mx, 3, 3, (1+POW_2(O/o_1))*rcc+2*(O/o_1)*rss ) ;

	/* Multiplies times -4 to obtain the Bloch matrix */
	gsl_matrix_scale ( red_mx, -4 ) ;

	return 0;
}		/* -----  end of function red_mat  ----- */


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

#pragma once

#ifndef texbatch_dfe6f3b4_430b_4424_98b6_d6b23235dc22_h
#define texbatch_dfe6f3b4_430b_4424_98b6_d6b23235dc22_h

#include <gslib/std.h>
#include <ariel/type.h>
#include <ariel/image.h>
#include <ariel/rectpack.h>
#include <ariel/rendersys.h>

__ariel_begin__

typedef render_texture2d texture2d;

#define _GS_BATCH_TEXTURE

#if defined(_GS_BATCH_TEXTURE)
#define location_key    texture2d*
#elif defined(_GS_BATCH_IMAGE)
#define location_key    const image*
#endif

class tex_batcher
{
public:
    typedef unordered_map<location_key, rectf> location_map;

public:
    tex_batcher();
    bool is_empty() const { return _rect_packer.is_empty(); }
    float get_width() const { return _rect_packer.get_width() + _gap; }
    float get_height() const { return _rect_packer.get_height() + _gap; }
    void add_image(const image* p);
    void add_texture(texture2d* p);
    void arrange();
    const location_map& get_location_map() const { return _location_map; }
    texture2d* create_texture(rendersys* rsys) const;
    void create_packed_image(image& img) const;
    void tracing() const;

protected:
    rect_packer         _rect_packer;
    location_map        _location_map;
    float               _gap;

private:
    void prepare_input_list(rp_input_list& inputs);
};

__ariel_end__

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib