text:
#ifndef SPARSE_AUC_AUC_OPT_METHODS_H
#define SPARSE_AUC_AUC_OPT_METHODS_H

#include <time.h>
#include <stdbool.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <limits.h>
// These are the third part library needed.
#include <cblas.h>
#include "fast_pcst.h"
#include "loss.h"

#define PI 3.14159265358979323846
#define sign(x) (x > 0) - (x < 0)
#define max(a, b) ((a) > (b) ? (a) : (b))
#define min(a, b) ((a) < (b) ? (a) : (b))
#define swap(a, b) { register double temp=(a);(a)=(b);(b)=temp; }
#define is_posi(x) ( x > 0.0 ? 1.0 : 0.0)
#define is_nega(x) ( x < 0.0 ? 1.0 : 0.0)

typedef struct {
    double *wt;
    double *wt_prev;
    double *aucs;
    double *rts;
    int auc_len; // how many auc evaluated.
    int total_iterations; // total iterations
    int total_epochs; // total epochs executed.
} AlgoResults;

typedef struct {
    int num_passes;
    int verbose;
    int step_len;
    int record_aucs;
    double stop_eps;
} GlobalParas;

AlgoResults *make_algo_results(int data_p, int total_num_eval);

bool free_algo_results(AlgoResults *re);

typedef struct {
    Array *re_nodes;
    Array *re_edges;
    double *prizes;
    double *costs;
    int num_pcst;
    double run_time;
    int num_iter;
} GraphStat;

typedef struct {
    const double *x_tr_vals;
    const int *x_tr_inds;
    const int *x_tr_poss;
    const int *x_tr_lens;
    const double *y_tr;
    bool is_sparse;
    int n;
    int p;
    // this is only for the graph operator.
    bool is_graph;
    int m; // number of edges.
    EdgePair *edges;
    double *weights;
    int g;
    double *proj_prizes;
    GraphStat *graph_stat;
} Data;

GraphStat *make_graph_stat(int p, int m);

bool free_graph_stat(GraphStat *graph_stat);

typedef struct {
    double val;
    int index;
} data_pair;

bool head_tail_binsearch(
        const EdgePair *edges, const double *costs, const double *prizes,
        int n, int m, int target_num_clusters, int root, int sparsity_low,
        int sparsity_high, int max_num_iter, PruningMethod pruning,
        int verbose, GraphStat *stat);

/**
 * SOLAM: Stochastic Online AUC Maximization
 * ---
 * BibTEX:
 * @inproceedings{ying2016stochastic,
 * title={Stochastic online AUC maximization},
 * author={Ying, Yiming and Wen, Longyin and Lyu, Siwei},
 * booktitle={Advances in neural information processing systems},
 * pages={451--459},
 * year={2016}
 * }
 * @param data
 * @param paras
 * @param re
 * @param para_xi
 * @param para_r
 * @author --- (Email: ---)
 * @return
 */
bool _algo_solam(Data *data,
                 GlobalParas *paras,
                 AlgoResults *re,
                 double para_xi,
                 double para_r);

/**
 * This function implements the algorithm proposed in the following paper.
 * Stochastic Proximal Algorithms for AUC Maximization.
 * ---
 * @inproceedings{natole2018stochastic,
 * title={Stochastic proximal algorithms for AUC maximization},
 * author={Natole, Michael and Ying, Yiming and Lyu, Siwei},
 * booktitle={International Conference on Machine Learning},
 * pages={3707--3716},
 * year={2018}}
 * ---
 * Do not use the function directly. Instead, call it by Python Wrapper.
 * @param data
 * @param paras
 * @param re
 * @param para_xi
 * @param para_l1_reg
 * @param para_l2_reg
 * @author --- (Email: ---)
 */
void _algo_spam(Data *data,
                GlobalParas *paras,
                AlgoResults *re,
                double para_xi,
                double para_l1_reg,
                double para_l2_reg);

/**
 * Stochastic Hard Thresholding for AUC maximization.
 * @param data
 * @param paras
 * @param re
 * @param para_s
 * @param para_b
 * @param para_c
 * @param para_l2_reg
 */
void _algo_sht_auc(Data *data,
                   GlobalParas *paras,
                   AlgoResults *re,
                   int version,
                   int operator_id,
                   int para_s,
                   int para_b,
                   double para_c,
                   double para_l2_reg);

void _algo_sto_iht(Data *data,
                   GlobalParas *paras,
                   AlgoResults *re,
                   int para_s,
                   int para_b,
                   double para_xi,
                   double para_l2_reg);

/**
 *
 * @param data
 * @param paras
 * @param re
 * @param para_s
 * @param para_tau
 * @param para_zeta
 * @param para_step_init
 * @param para_l2
 */
void _algo_hsg_ht(Data *data,
                  GlobalParas *paras,
                  AlgoResults *re,
                  int para_s,
                  double para_tau,
                  double para_zeta,
                  double para_step_init,
                  double para_l2);

int _hard_thresholding(double *arr, int n, int k);

/**
 *
 * @param data
 * @param paras
 * @param re
 * @param para_tau
 * @param para_eta
 * @param para_lambda
 */
void _algo_opauc(Data *data,
                 GlobalParas *paras,
                 AlgoResults *re,
                 int para_tau,
                 double para_eta,
                 double para_lambda);

/**
 * This function implements the algorithm, FSAUC, proposed in the following paper:
 * ---
 * @inproceedings{liu2018fast,
 * title={Fast stochastic AUC maximization with O (1/n)-convergence rate},
 * author={Liu, Mingrui and Zhang, Xiaoxuan and Chen, Zaiyi and Wang, Xiaoyu and Yang, Tianbao},
 * booktitle={International Conference on Machine Learning},
 * pages={3195--3203},
 * year={2018}}
 * ---
 * @param data
 * @param paras
 * @param re
 * @param para_r
 * @param para_g
 */
void _algo_fsauc(Data *data,
                 GlobalParas *paras,
                 AlgoResults *re,
                 double para_r,
                 double para_g);

#endif //SPARSE_AUC_AUC_OPT_METHODS_H

max_stars_repo_name: baojianzhou/sparse-auc

max_issues_repo_name: baojianzhou/sparse-auc

max_forks_repo_name: baojianzhou/sparse-auc


text:
//gcc -g meanpreds.c -o meanpreds
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "my_misc.c"
#define isnum(c) (((c>='0') && (c<='9'))||(c=='-')||(c=='+')||(c==' '))
//static char pr[1024];
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix_double.h>
//#include <gsl/gsl_multimin.h>
#include <gsl/gsl_multiroots.h>
//#include "my_misc.c"
//#define GSL 1
#undef GSL
#define GSL 1
//#define GSL 0
int GSLtag;
void least_gsl_error_handler( const char * reason,
			      const char * file,
			      int line,
			      int gsl_errno)
{
  //  printf("%s:%d: %s (gsl_error: %s)\n",
  //	 file, line, reason, gsl_strerror(gsl_errno));
  printf("%s:%d: %s (gsl_error %d)\n",
	 file, line, reason,gsl_errno);
  GSLtag=1;
}
double VarY;
void calc_Ainvb(double *M, double *A_data, double *b_data, int nx, int ndata,int nxmarume)
{
  /* b= A M  --> M= Ainv b
     b  in R^{ndata x 1}
     A  in R^{ndata x nx}   (ndata>nx)
     M  in R^{nx    x 1}
   */
  int i;
  gsl_vector *S;
  gsl_vector *work;
  gsl_vector *x;
  gsl_vector *b;
  gsl_matrix *V;
  gsl_matrix *A;
  if(ndata<=0){
    M[0]=1;
    for(i=1;i<nx;i++) M[i]=1./nx;
  }
  else{
    A =gsl_matrix_alloc(ndata,nx);
    b =gsl_vector_alloc(ndata);b->data=b_data;
    A->data=A_data;
    //  for(i=0;i<ndata;i++){printf("%e %e %d #bi,Ai\n",A->data[i],b->data[i],i); }

    V = gsl_matrix_alloc (nx,nx);
    S = gsl_vector_alloc (nx);
    work = gsl_vector_alloc(nx);
    x = gsl_vector_alloc(nx);

    gsl_set_error_handler( &least_gsl_error_handler );GSLtag=0;

#if GSL == 1
    gsl_linalg_SV_decomp_jacobi (A, V, S); //higer accuracy than Golub-Reinsh
#elif GSL == 2
    gsl_linalg_SV_decomp (A, V, S, work); //Golub-Reinsh
#elif GSL == 3
    {//faster
      gsl_matrix *X;
      X = gsl_matrix_alloc (nx,nx);
      gsl_linalg_SV_decomp_mod (A, X, V, S, work);
      gsl_matrix_free(X);
    }
#endif
    //    for(i=nxmarume;i<nx;i++){

    fprintf(stdout,"Lambda");
    for(i=0;i<nx;i++){
      //      if(nx>40) break;
      //      fprintf(stdout,"S->data[%d]=%e;",i,S->data[i]);
      fprintf(stdout,"(%d)%e",i,S->data[i]);
    }
    fprintf(stdout,"\n");
    VarY=0;
    if(GSLtag==1){
      fprintf(stderr,"\nEigenValue=");for(i=0;i<nx;i++) fprintf(stderr,"%e ",S->data[i]); fprintf(stderr,"\n");
      for(i=0;i<nx;i++) M[i]=1./nx;
    }
    else{
      if(nxmarume>=0){//Principal Component Analysis
	for(i=nxmarume;i<nx;i++) S->data[i]=0;
	for(i=0;i<nxmarume;i++) VarY+=S->data[i];
	gsl_linalg_SV_solve(A,V,S,b,x);
	for(i=0;i<nx;i++) M[i]=gsl_vector_get(x,i);
	//    {//check
	//      double err=0,yhat;
	//      int j;
	//      for(i=0;i<ndata;i++){
	//	if(i>10) break;
	//	yhat=b_data[i];
	//	for(j=0;j<nx;j++){
	//	  yhat-= M[j]*A_data[i*nx+j];
	//	}
	//	err+=square(yhat);
	//	fprintf(stderr,"err=%e=%e/%d. yhat=%e=(%e",err/(i+1.),err,i,yhat,b_data[i]);
	//	for(j=0;j<nx;j++) fprintf(stderr,"-(%e)*(%e)",M[j],A_data[i*nx+j]);
	//	fprintf(stderr,")^2\n");
	//      }
	//    }
      }
      else{//Minor Component Analysis
	int ii=0;
	int im=-nxmarume;
	for(i=nx-1;i>=0;i--){
	  if(S->data[i]>1e-20){
	    ii++;
	    VarY+=S->data[i];
	    if(ii>=im) break;
	  }
	}
	int j;

	double M0=-V->data[i];
	for(j=1;j<nx;j++){
	  M[j-1]=V->data[j*nx+i]/M0;
	}
//	double M0=-V->data[i*nx];
//	for(j=1;j<nx;j++){
//	  M[j-1]=V->data[i*nx+j]/M0;
//	}

      }
    }
    gsl_vector_free(S);
    gsl_vector_free(work);
    gsl_vector_free(x);
    gsl_matrix_free(V);
    gsl_matrix_free(A);
    gsl_vector_free(b);
  }
}
void calc_AtinvbMCA(double *M, double *At_data, double *b_data, int nx, int ndata,int nxmarume)
{
  /* b= A M  --> M= Ainv b
     b  in R^{ndata x 1}
     A=At  in R^{ndata x nx}
     M  in R^{nx    x 1}
   */
  int nx1=nx+1,i,j,jnx1;
  gsl_matrix *A  =gsl_matrix_alloc(ndata,nx1);
  for(j=0;j<ndata;j++){
    jnx1=j*nx1;
    A->data[jnx1]=b_data[j];
    for(i=0;i<nx;i++){
      A->data[jnx1+i+1]=At_data[i*ndata+j];
    }
  }
  calc_Ainvb(M, (double *)(A->data), b_data, nx1, ndata, -nxmarume);

  gsl_matrix_free(A);
}
void calc_Atinvb(double *M, double *At_data, double *b_data, int nx, int ndata,int nxmarume)
{
  /* b= A M  --> M= Ainv b
     b  in R^{ndata x 1}
     A=At  in R^{ndata x nx}
     M  in R^{nx    x 1}
   */
  gsl_matrix *At =gsl_matrix_alloc(nx,ndata);
  gsl_matrix *A  =gsl_matrix_alloc(ndata,nx);
  At->data=At_data;
  gsl_matrix_transpose_memcpy(A, At);
  calc_Ainvb(M, (double *)(A->data), b_data, nx, ndata, nxmarume);
  gsl_matrix_free(A);
  gsl_matrix_free(At);
}
void calc_yp(double *y, double *a_data, int nx, int ndata,int nxmarume)
{//M[i]=Ainv b
  int i, l,j;
  double sumS=0,sumSj=0;
  gsl_vector *S;
  gsl_vector *work;
  gsl_vector *x;
  //  gsl_vector *b;
  gsl_matrix *V;
  gsl_matrix *A;
  A =gsl_matrix_alloc(ndata,nx);
  A->data=a_data;
  //    b =gsl_vector_alloc(ndata);
  //    b->data=b_data;
  V = gsl_matrix_alloc (nx,nx);
  S = gsl_vector_alloc (nx);
  work = gsl_vector_alloc(nx);
  x = gsl_vector_alloc(nx);

  gsl_set_error_handler( &least_gsl_error_handler );GSLtag=0;

#if GSL == 1
  gsl_linalg_SV_decomp_jacobi (A, V, S); //higer accuracy than Golub-Reinsh
#elif GSL == 2
  gsl_linalg_SV_decomp (A, V, S, work); //Golub-Reinsh
#elif GSL == 3
  {//faster
    gsl_matrix *X;
    X = gsl_matrix_alloc (nx,nx);
    gsl_linalg_SV_decomp_mod (A, X, V, S, work);
    gsl_matrix_free(X);
  }
#endif
  //    for(i=nxmarume;i<nx;i++){
  sumS=sumSj=0;
  for(l=0;l<nx;l++) sumS+=gsl_vector_get(S,l);
  if(nxmarume>=0){
    for(i=0;i<ndata;i++){
      y[i]=0;
      for(l=0;l<nxmarume;l++){
	//	fprintf(stderr,"S%d=%e.\n",l,gsl_vector_get(S,l));
	for(j=0;j<nx;j++){
	  //	y[i]+=S->data[l]*A->data[i*nx+l]*V->data[j*nx+l];
	  y[i]+=gsl_vector_get(S,l)*gsl_matrix_get(A,i,l)*gsl_matrix_get(V,j,l);
	}
      }
      y[i]/=nx;
      if(i<5){
	sumSj+=gsl_vector_get(S,i);
	fprintf(stderr,"[%d]y=%e, s=%e, CumulativeProp=%f=%f/%f\n",i,y[i],gsl_vector_get(S,i),sumSj/sumS,sumSj,sumS);
      }
    }
  }
  else{
    int ll=0;
    for(i=0;i<ndata;i++){
      y[i]=0;
      for(l=nx;;l--){
	if(gsl_vector_get(S,l)<=1e-30) continue;
	if(i==0) fprintf(stderr,"S%d=%e.\n",l,gsl_vector_get(S,l));
	for(j=0;j<nx;j++){
	  y[i]+=gsl_vector_get(S,l)*gsl_matrix_get(A,i,l)*gsl_matrix_get(V,j,l);
	}
	if(++ll>=-nxmarume) break;
      }
      y[i]/=nx;
      if(i<5){
	sumSj+=gsl_vector_get(S,i);
	if(i<5) fprintf(stderr,"[%d]y=%e, s=%e, CumulativeProp=%f=%f/%f\n",i,y[i],gsl_vector_get(S,i),sumSj/sumS,sumSj,sumS);
      }
    }
  }
  gsl_vector_free(S);
  gsl_vector_free(work);
  gsl_vector_free(x);
  gsl_matrix_free(V);
}

////For Maximizing the Entropy of the training data
////For Minimizing the minus of the Entropy of the training data
//Lagrange Method
int
my_df_Entropy(const gsl_vector *v, void *params, gsl_vector *df)
{
  double **dp = (double **)params;
  int nd=(int)((double)*dp[0]);
  int nl=(int)((double)*dp[1]);
  double *y=(double *)dp[2];
  double *f=(double *)dp[3];
  int sgn=(int)((double)*dp[4]);
  double *yp=(double*)malloc(sizeof(double)*nd);
  double *e2=(double*)malloc(sizeof(double)*nd);
  double *e =(double*)malloc(sizeof(double)*nd);
  double *dHdp=(double*)malloc(sizeof(double)*nd);
  double *dpde2=(double*)malloc(sizeof(double)*nd*nd);
  double E2,lognd=log(nd),pi;
  int i,j,k,l;

  E2=0;
  for(i=0;i<nd;i++){
    yp[i]=0; for(j=0;j<nl;j++) yp[i]+=square(v->data[j])*f[j*nd+i];
    e[i]=(y[i]-yp[i]);
    e2[i]=e[i]*e[i];
    E2+=e2[i];
  }
  for(i=0;i<nd;i++){
    pi=e2[i]/E2;
    dHdp[i]=(log(pi)+1.)/lognd;//for the minus Entropy
    for(k=0;k<nd;k++){
      if(i == k) dpde2[i*nd+k]=1./E2-e2[i]/E2/E2;
      else dpde2[i*nd+k]=-e2[i]/E2/E2;
    }
  }
  //  double dfnorm=0;
  //  fprintf(stdout,"df=");
  for(l=0;l<nl;l++){
    df->data[l]=0;
    for(i=0;i<nd;i++){
      for(k=0;k<nd;k++){
	df->data[l]+=dHdp[i]*dpde2[i*nd+k]*(2.*e[k]*(-2.*v->data[l]*f[l*nd+k]));
	//	df->data[l]+=dHdp[i]*dpde2[i*nd+k]*(-2.*e[k]*f[l*nd+k]);
      }
      df->data[l]+= 2.*v->data[nl]*v->data[l]*sgn;//for Penalty
      //      df->data[l]-= 2.*v->data[nl]*v->data[l];//for Penalty
      //      df->data[l] = C0*(-2.*e[i]*f[l*nd+i])/E0 +C1*df->data[l];//for H=C0*E2/E0+C1*Entropy
    }
    //    dfnorm+=square(df->data[l]);
    //    fprintf(stdout,"%.4f ",df->data[l]);
  }
  //Lagrange for sumv=1. v->data[nd]=lambda
  double sumv=0;//for Lagrange
  for(j=0;j<nl;j++) sumv+=square(v->data[j]);//for Lagrange
  df->data[nl]=(sumv-1.)*sgn;//for Lagrange L=H+lambda*(sum c_l-1)
  //  df->data[nl]=1.-sumv;//for Lagrange L=H-lambda*(sum c_l-1)
  //  fprintf(stdout,";dfnorm=%.4e\n",dfnorm);
  free(yp);
  free(e);
  free(e2);
  return GSL_SUCCESS;
}
int
gsl_print_state (size_t iter, gsl_multiroot_fsolver * s, void *params)
{
  double **dp = (double **)params;
  int nd=(int)((double)*dp[0]);
  int nl=(int)((double)*dp[1]);
  double *y=(double *)dp[2];
  double *f=(double *)dp[3];
  double *yp=(double*)malloc(sizeof(double)*nd);
  double *e2=(double*)malloc(sizeof(double)*nd);
  double *e =(double*)malloc(sizeof(double)*nd);
  double E2,pi;
  int i,j;

  E2=0;
  for(i=0;i<nd;i++){
    yp[i]=0; for(j=0;j<nl;j++) yp[i]+=square(gsl_vector_get (s->x, j))*f[j*nd+i];
    e[i]=(y[i]-yp[i]);
    e2[i]=e[i]*e[i];
    E2+=e2[i];
  }
  double H=0;//
  for(i=0;i<nd;i++){
    pi=e2[i]/E2;
    H += pi*log(pi);//the minus Entropy
  }
  H/=log(nd);
  printf ("iter = %3u H=%f x0=%.3e...%.3e "
	  "dL/dx=%.3e...%3e\n",
	  (unsigned int)iter,-H,
	  square(gsl_vector_get (s->x, 0)),
	  gsl_vector_get (s->x, nl),
	  gsl_vector_get (s->f, 0),
	  gsl_vector_get (s->f, nl));
  return 0;
}

//////For Maximizing the Entropy of the training data
//////For Minimizing the minus of the Entropy of the training data
//double
//my_f_Entropy(const gsl_vector *v, void *params)
//{
//  //params
//  /*
//    e_i=y_i- sum_l v[l] f[l][i] (l=0:nx-1)
//    p_i= e_i^2/ sum_j e_j^2
//    H= sum_i p_i log(p_i);
//
//  */
//  double **dp = (double **)params;
//  int nd=(int)((double)*dp[0]);
//  int nl=(int)((double)*dp[1]);
//  double *y=(double *)dp[2];
//  double *f=(double *)dp[3];
//  double C0=(double)*dp[4];
////  double C1=(double)*dp[5];
////  double E0=(double)*dp[6];
//  double *e2=(double*)malloc(sizeof(double)*nd);
//  double E2;
//  double H,ypi,pi;
//  int i,j;
//  E2=0;
//  for(i=0;i<nd;i++){
//    ypi=0; for(j=0;j<nl;j++) ypi+=square(v->data[j])*f[j*nd+i];
//    e2[i]=square(ypi-y[i]);
//    E2+=e2[i];
//  }
//  H=0;
//  for(i=0;i<nd;i++){
//    pi=e2[i]/E2;
//    H += pi*log(pi);//the minus Entropy
//  }
//  H/=log(nd);
//
//  //  H=C0*E2/E0+C1*H;//for H=C0*E2/E0+C1*Entropy
//
//  //Lagrange for sumv=1. v->data[nd]=lambda
//  double sumv=0;//for Lagrange
//  for(j=0;j<nl;j++) sumv+=square(v->data[j]);//for Lagrange
//  //  if(C0*(sumv-1.)<0) C0*=-1;
//  //  H+=v->data[nl]*(sumv-1.);//for Lagrange
//  fprintf(stderr,"H=%f,sumv=%e",H,sumv);
//  H+=C0*(sumv-1.);//for Penalty method
//  /////
//  free(e2);
//  return H;
//}
//
///* The gradient of f, df = (df/dx, df/dy). */
//void
//my_df_Entropy(const gsl_vector *v, void *params,
//       gsl_vector *df)
//{
//  double **dp = (double **)params;
//  int nd=(int)((double)*dp[0]);
//  int nl=(int)((double)*dp[1]);
//  double *y=(double *)dp[2];
//  double *f=(double *)dp[3];
//  double C0=(double)*dp[4];
//  //  double C1=(double)*dp[5];
//  //  double E0=(double)*dp[6];
//  double *yp=(double*)malloc(sizeof(double)*nd);
//  double *e2=(double*)malloc(sizeof(double)*nd);
//  double *e =(double*)malloc(sizeof(double)*nd);
//  double *dHdp=(double*)malloc(sizeof(double)*nd);
//  double *dpde2=(double*)malloc(sizeof(double)*nd*nd);
//  double E2,lognd=log(nd),pi;
//  int i,j,k,l;
////  {
////    double sumv;
////    for(j=0;j<nl;j++) sumv+=square(v->data[j]);//for Lagrange
////    if(C0*(sumv-1.)<0) C0*=-1;
////  }
//  E2=0;
//  for(i=0;i<nd;i++){
//    yp[i]=0; for(j=0;j<nl;j++) yp[i]+=square(v->data[j])*f[j*nd+i];
//    e[i]=(y[i]-yp[i]);
//    e2[i]=e[i]*e[i];
//    E2+=e2[i];
//  }
//  for(i=0;i<nd;i++){
//    pi=e2[i]/E2;
//    //    dHdp[i]=(-log(pi)-1.)/lognd;
//    dHdp[i]=(log(pi)+1.)/lognd;//for the minus Entropy
//    for(k=0;k<nd;k++){
//      if(i == k) dpde2[i*nd+k]=1./E2-e2[i]/E2/E2;
//      else dpde2[i*nd+k]=-e2[i]/E2/E2;
//    }
//  }
//  double dfnorm=0;
//  //  fprintf(stdout,"df=");
//  for(l=0;l<nl;l++){
//    df->data[l]=0;
//    for(i=0;i<nd;i++){
//      for(k=0;k<nd;k++){
//	df->data[l]+=dHdp[i]*dpde2[i*nd+k]*(2.*e[k]*(-2.*v->data[l]*f[l*nd+k]));
//	//	df->data[l]+=dHdp[i]*dpde2[i*nd+k]*(-2.*e[k]*f[l*nd+k]);
//      }
//      df->data[l]+= 2.*C0*v->data[l];//for Penalty
//      //      df->data[l] = C0*(-2.*e[i]*f[l*nd+i])/E0 +C1*df->data[l];//for H=C0*E2/E0+C1*Entropy
//    }
//    dfnorm+=square(df->data[l]);
//    //    fprintf(stdout,"%.4f ",df->data[l]);
//  }
//  //Lagrange for sumv=1. v->data[nd]=lambda
//  //  double sumv=0;//for Lagrange
//  //  for(j=0;j<nl;j++) sumv+=square(v->data[j]);//for Lagrange
//  //  df->data[nl]=sumv-1.;//for Lagrange
//
//  fprintf(stdout,";dfnorm=%.4e\n",dfnorm);
//  free(yp);
//  free(e);
//  free(e2);
//}
//
///* Compute both f and df together. */
//void
//my_fdf_Entropy(const gsl_vector *x, void *params,
//	double *f, gsl_vector *df)
//{
//  *f = my_f_Entropy(x, params);
//  my_df_Entropy(x, params, df);
//}

////
typedef struct{
  int r1;
  int r2;
  double r3;
  int nr;
  double r12;
  double *r;
  double ymin;
  double ymax;
} RESOLUTION;
#define ZERO 1e-20
int init_res(RESOLUTION *res)
{
  int i;
  char buff[32];
  if(res->r2==0) return(-1);
  if(fabs(res->ymax-res->ymin)<ZERO) {
    res->r2=0;
    return(-1);
  }
  if(res->r3<-ZERO){//negrect resolution for the meanpred
    res->r2=0;
    return(-1);
  }
  res->r12=(double)res->r1/res->r2;
  //  res->nr=(res->ymax-res->ymin)/res->r12;
  res->nr=(res->ymax-res->ymin)/res->r12+0.5;//??
  //  res->r12=(res->ymax-res->ymin)/res->nr;//for ymin and ymax more reliable than r1 and r2 for ijcnn06
  res->r=(double*)malloc(sizeof(double)*(res->nr+1));
  if(fabs(res->r3)>ZERO){
    for(i=0;i<=res->nr;i++){
      sprintf(buff,"%.10le",((int)((i*res->r12+res->ymin)/res->r3+0.5))*res->r3);//r3$B7e$G;M<N8^F~(B
      sscanf(buff,"%lf",&res->r[i]);
    }
  }
  else{
    for(i=0;i<=res->nr;i++){
      sprintf(buff,"%.10le",(double)i*res->r12+res->ymin);
      sscanf(buff,"%lf",&res->r[i]);
    }
  }
  return(0);
}
double y_res(double y,RESOLUTION *res)
{
  int ii;
  if(res->r2==0) return(y);
  //  ii=(int)((y-ymin)/r12+0.5);//4sha5nyu
  ii=floor((y-res->ymin)/res->r12+0.5);//4sha5nyu
  if(y<res->ymin) return((double)res->r[0]);
  if(y>res->ymax) return((double)res->r[res->nr]);
  return((double)res->r[ii]);
}
#define buffsize 1024
int main(int argc, char **argv)
{
  int i,j,l;
  FILE *fp;
  char buff[buffsize];
  double *yp,_yp;
  double *y,_y,*yr,_yr;
  int *n,_n;
  int num=0,num0=0;
  //  char *fnpredict="predict+.dat";
  //  char *fnloss="loss+.dat";
  int intmode=0, j0=1,nfiles=argc-1;
  char fn_is[256],fn_pred[256],fn_pred0[256],fn_prob[256];
  double _is,_mseis,_n_msetrain,_msetrain,_n_cells2,_nCells1,_n_mse,_mse,_msetrainis;
  double *prob,*sigma2,meansigma2,probtotal;
  FILE *fp0;
  int bayesmode=0; //0 (for flat probability) or 1 (for probatility decaying according to MSE).
  double bayesopt1=0;//for adjusting probability (1 for same probatility)
  double bayesopt2=0;//0 for original L, 1 for L_derivative requiring nfiles>=3.
  double bayesopt3=0;//??
  double bayesopt4=0;//??
  //  double *ds;
  //  int *s;
  double Lvar,Lvarval,Lval,LD=0,_LD;//
  double skew,kurtosis,entropy,skew2,skewa;//skew=3rd moment, kurtosis=4th moment
  double m2,m3,m4,m5,m6,m3a,m5a;//5th moment, 6th moment
  double m3p,m3m;//m3plus,m3minus
  int nm3p,nm3m;
  double mcov;//mean covariance
  double Lvar0;
  double Ltrain,LtrainE,H0;
  double *var,*vv;
  int *cc,ccmax=0;
  double _var,_vv,_cc;
  double *kyorix2,_kyorix2,_e2;
  int LDmode=2;
  RESOLUTION res={0,0,0,};
  double *yp00;
  int nop=0;
  if(argc<2){
    fprintf(stderr,"Usage: %s [ib:<intmode>:<basesmode>[:<probfile>]] <pred&is-body#1> <pred&is-body#1> ... \n",argv[0]);
    fprintf(stderr,"<intmode>:1 for intmode, 0 for real.\n");
    fprintf(stderr,"<bayesmode>:0 simple ensemble.?\n");
    fprintf(stderr,"<bayesmode>:1 for bayes mode1 exp.\n");
    fprintf(stderr,"<bayesmode>:2 for bayes mode2 YpInvY.\n");
    fprintf(stderr,"<bayesmode>:-1 for reading from probfile.\n");
    fprintf(stderr,"res:<r1>:<r2>:<ymin>:<ymax>.\n");
    fprintf(stderr,"Old Usage: %s [<int>:<bases>] <pred&is-body#1> <pred&is-body#1> ... \n",argv[0]);
    fprintf(stderr,"par:paramfile.\n");
    exit(1);
  }
  //  fprintf(stderr,"last arg[%d]='%s'\n",argc,argv[argc-1]);
  int parfile=0;
  //
  for(i=1;i<argc;i++){
    if(strncmp(argv[i],"par:",4)==0){//for compatibility to old version
      parfile=1;
      break;
    }
  }
  char *_argv0;
  if(parfile==1){
    char *fn=&argv[i][4];
    FILE *fp;
    if((fp=fopen(fn,"r"))==NULL){
      fprintf(stderr,"Param file (%s) open error.\n",fn);
      return(-1);
    }

    fseek( fp, 0L, SEEK_END );
    int fsize=ftell(fp);
    fseek( fp, 0L, SEEK_SET );//rewind(fp);//
    //    fprintf(stderr,"fsize=%d\n",fsize);
    _argv0=(char*)malloc(sizeof(char)*fsize);
    fgets(_argv0,fsize,fp);
    fclose(fp);
    argc=1;
    char *p=_argv0;
    for(j=0;j<fsize;j++) if(p[j]==' ') argc++;
    argv=(char**)malloc(sizeof(char*)*(argc));
    for(j=1;j<argc;j++){
      for(p++;;p++) if(*p==' ') {*p=0;break;}
      for(p++;;p++) if(*p!=' ') break;
      argv[j]=p;
    }
  }
  intmode=0;
  bayesmode=0;
  j0=1;
  nfiles=argc-1;
  for(i=j0;i<argc;i++){
    if(strncmp(argv[i],"int",3)==0){//for compatibility to old version
      intmode=1;
      j0++;
      nfiles--;
    }
    else if(strncmp(argv[i],"ib:",3)==0){
      sscanf(&(argv[i][3]),"%d:%d:%lf:%lf:%lf:%lf",&intmode,&bayesmode,&bayesopt1,&bayesopt2,&bayesopt3,&bayesopt4);
      if(bayesmode==-1) sprintf(fn_prob,"%s",&(argv[1][8]));
      j0++;
      nfiles--;
    }
    else if(strncmp(argv[i],"ry:",3)==0){
      sscanf(&argv[i][3],"%d:%d:%lf:%lf:%lf",&res.r1,&res.r2,&res.r3,&res.ymin,&res.ymax);
      j0++;
      nfiles--;
      init_res(&res);
    }
    else if(strncmp(argv[i],"LDm:",4)==0){
      sscanf(&argv[i][4],"%d",&LDmode);
      j0++;
      nfiles--;
      init_res(&res);
    }
    else if(strncmp(argv[i],"nop:",4)==0){
      sscanf(&argv[i][4],"%d",&nop);
      j0++;
      nfiles--;
    }
  }
  ////
  if(nop!=1){
    printf("#start:%s",argv[0]);
    for(i=1;i<argc;i++){
      printf(" %s",argv[i]);
    }
  }
//  else{//if(nop==1){
//    fprintf(stdout,"+");
//    fflush(stdout);
//  }

  int (*net_printf)(char *);
  if(nop==0) net_printf=printf1;
  else net_printf=noprintf;

  prob =(double*) malloc(sizeof(double)*nfiles);
  sigma2=(double*) malloc(sizeof(double)*nfiles);
  meansigma2=0;
  for(j=0;j<nfiles;j++) argv[j+j0][strlen(argv[j+j0])-8]=0; //argv[j+j0]=fn_bodypred.dat
  for(j=0;j<nfiles;j++){
    sprintf(fn_is,"%sis.dat",argv[j+j0]);
    if((fp=fopen(fn_is,"r"))==NULL){
      fprintf(stderr,"#Execute without '%s'\n",fn_is);
      //      fprintf(stderr,"file(%s) could not be opened!\n",fn_is);
      num=0; sigma2[j]=0; meansigma2+=0;//050914
      num0=0;
    }
    else{
      fscanf(fp,"%lf%lf%lf%lf%lf%lf%lf%lf%lf",&_is,&_mseis,&_n_msetrain,&_msetrain,&_n_cells2,&_nCells1,&_n_mse,&_mse,&_msetrainis);
      //      num=_n_msetrain; sigma2[j]=_msetrain; meansigma2+=_msetrain;//050902??
      //      num=_n_mse; sigma2[j]=_mse; meansigma2+=_mse;//050823
      num=_n_mse; sigma2[j]=_msetrain; meansigma2+=_msetrain;//050914
      num0=_n_msetrain;
      //      fprintf(stderr,"sigma2[%d]%e\n",j,sigma2[j]);
      fclose(fp);
    }
  }
  meansigma2/=nfiles;
  ////////////////////////////////
  if(bayesmode>=100){}//end of bayesmode>=100
  else{//bayesmode <100
    if(bayesmode==0){//same probability
      probtotal=nfiles;
      for(j=0;j<nfiles;j++)  prob[j]=1./probtotal;
      probtotal=0; for(j=0;j<nfiles;j++) probtotal+=prob[j];
    }
    else if(bayesmode==1){//probatility depending on MSE(=sigma2)
      int jj;double probj;//probj is 1/prob[j] for avoiding marumegosa.
      probtotal=0;
      for(j=0;j<nfiles;j++){
	probj=0;
	for(jj=0;jj<nfiles;jj++){
	  if(jj==j) probj+=1;
	  else probj+= exp((sigma2[j]-sigma2[jj])*bayesopt1/meansigma2);
	}
	prob[j]=1./probj;
	probtotal+=prob[j];
      }
      for(j=0;j<nfiles;j++)  prob[j]/=probtotal;
      probtotal=0; for(j=0;j<nfiles;j++) probtotal+=prob[j];
    }
    ////////////////////////////////
    int m_LD=bayesopt2,nfiles2;//bayesopt2=m for LD= Loss_D = sum_{j=1}^m ((f_{l0+j}-f_{l0})^2 +(f_{l0-j}-f_{l0})^2)
    if(bayesmode==3) m_LD=0;
    else m_LD=bayesopt2;
    if(nfiles<m_LD*2+1) m_LD=(nfiles-1)/2;// (nfiles,m_LD)=(2,0),(3,1),(4,1),(5,2),(6,2);
    nfiles2=nfiles-m_LD*2;//for LD of generalization error
    int nmean=nfiles-nfiles2+1;//nmean=2*m_LD+1
    //    double *ww =(double*) malloc(sizeof(double)*nfiles*(nmean));
    //#define LDmode 2
    //#define LDmode 4
    //	//LDmode==4 is only for ib:0:0:0:m ???
    //#if LDmode == 2
    //    double *ww =(double*) malloc(sizeof(double)*(nfiles+1)*(nmean));
    //#elif  LDmode == 4
    //    double *ww0 =(double*) malloc(sizeof(double)*(nfiles+1)*(nmean)*2);
    //    double *ww =(double*) malloc(sizeof(double)*(nfiles+1)*(nmean)*(nmean));
    double *ww =(double*) malloc(sizeof(double)*(nfiles)*(nmean)*(nmean));
    //#endif
    //    double *ww=&ww0[1];
    double *ww0=&ww[(m_LD+m_LD*nmean)*nfiles];

    {//////////for Lvar0
      //      double *yp0=(double*) malloc(sizeof(double)*num0*nfiles);
      //      double *yt0=(double*) malloc(sizeof(double)*(nfiles+1)*num0+2);
      //      double *yt=(double*) &yt0[2];//yt[0]=num0,yt[1]=nmean;

      double *yt=(double*) malloc(sizeof(double)*(nfiles+1)*num0);
      double *yp01=(double*) &yt[num0];
      double *yp =(double*) malloc(sizeof(double)*num0);
      //      nfiles2=nfiles-2*bayesopt2;//for LD of generalization error
      if(nfiles<3) {m_LD=bayesopt2=0;nfiles2=nfiles;}//training ?
      //      else {m_LD=1;nfiles2=nfiles-2*bayesopt2;}//for LD of generalization error
      //      probtotal=0;
      //      for(j=m_LD;j<nfiles2+m_LD;j++){//prediction of training data?
      for(j=0;j<nfiles;j++){//prediction of training data?
#define NEWp0
#ifdef NEWp0
	{
	  char *p1,*p;
	  sprintf(fn_pred0,"%s",argv[j+j0]);
	  for(p=fn_pred0;;p++){if(*p=='+') break;}
	  p1=++p;
	  for(;;p++){if(*p=='+') break;}
	  ++p;
	  sprintf(p1,"%spred0.dat",p);
	}
#else
	sprintf(fn_pred0,"%spred0.dat",argv[j+j0]);
#endif
	fp0=fopen(fn_pred0,"r");
	for(i=0;i<num0;i++){
	  fgets(buff,buffsize,fp0);
	  sscanf(buff,"%lf%d%lf%lf%lf%lf%lf",&_yp,&_n,&_y,&_yr,&_var,&_cc,&_vv);//sscanf(buff,"%lf%d%lf%lf",&_yp,&_n,&_y,&_yr);
	  if(res.r3<0) _yp=_yr; //neglect resolution for the meanpred
	  else _yp=y_res(_yp,&res);
	  //	  yp0[i*nfiles+j] =_yp;
	  yp01[j*num0+i] =_yp;
	  yt[i]=_y;
	}
	fclose(fp0);
      }
      if(bayesmode==2){
	//	if(1==0 && nfiles==1) ww[0]=1;
	//	else
	{
	  for(j=0;j<nmean;j++){
	    for(l=0;l<nmean;l++){
	      //LDmode = 2,4,8 neighbours
	      if(LDmode==2 && l!=m_LD) continue;//only yoko
	      if(LDmode==4 && (l-m_LD)*(j-m_LD)!=0 ) continue;//only tate&yoko
	      //calc_AtinvbMCA((double *)(&ww[j*nfiles]),(double *)(&yp01[j*num0]),(double *)(&yt[0]),(int)(nfiles2),(int)num0,(int)bayesopt1);
	      calc_Atinvb((double *)(&ww[(l*nmean+j)*nfiles]),(double *)(&yp01[j*num0]),(double *)(&yt[0]),(int)(nfiles2+l-m_LD),(int)num0,(int)bayesopt1);
	      //	  calc_Ainvb((double *)(&ww[j*nfiles]),(double *)(&yp0[j]),(double *)(&yt[j]),(int)(nfiles2),(int)num0,(int)bayesopt1);
	      //	  calc_Ainvb((double *)(&ww[j*nfiles]),(double *)(&yp01[j*num0]),(double *)(&yt[0]),(int)(nfiles2),(int)num0,(int)bayesopt1);
	      //	    calc_Ainvb((double *)(&ww[j*nfiles]),(double *)(&yp01[j*num0]),(double *)(&yt[0]),(int)(nfiles2),(int)num0,(int)bayesopt1);
	      //	    for(i=0;i<=nfiles2;i++) ww[j*nfiles+i]=1./nfiles2;
	      //	    for(i=0;i<10;i++) fprintf(stderr,"j%2d,i%2d)yt%+e yp%+e %+e\n",j,i,yt[i],yp01[(j)*num0+i],yp01[(j+1)*num0+i]);
	      if(fabs(VarY)<1e-10){
		for(i=0;i<nfiles2+l-m_LD;i++) ww[(l*nmean+j)*nfiles+i]=1./(nfiles2+l-m_LD);
	      }
	      fprintf(stderr,"w[%d,%d]:",j,l);
	      for(i=0;i<nfiles2+l-m_LD;i++) fprintf(stderr,"%+.4f ",ww[(l*nmean+j)*nfiles+i]);
	      fprintf(stderr,"\n");
	    }
	  }
	}
      }

      else if(bayesmode==3){//ib:0:3:10:0:0 ib:0:3:itermax:0:0 H=C0*E2/E0+C1*Entropy
	for(j=0;j<nmean;j++){
	  for(l=0;l<nmean;l++){
	    if(LDmode==2 && l!=m_LD) continue;//only yoko
	    if(LDmode==4 && (l-m_LD)*(j-m_LD)!=0 ) continue;//only tate&yoko
	    {
	      size_t iter = 0;
	      size_t itermax=bayesopt1;if(itermax<1e-20) itermax=100;
	      int status;
	      double nl=(double)(nfiles2+l-m_LD);
	      double nd=(double)num0;
	      int nl1=nl+1;
	      double *par[7];
	      double sgn;
	      if(bayesopt2>=0) sgn=1; else sgn=-1;
	      par[0]=(double *)&nd;
	      par[1]=&nl;
	      par[2]=&yt[0];
	      par[3]=&yp01[j*num0];
	      par[4]=&sgn;
	      const gsl_multiroot_fsolver_type *T;
	      gsl_multiroot_fsolver *s;
	      gsl_multiroot_function df = {&my_df_Entropy, nl1, &par};
	      gsl_vector *x = gsl_vector_alloc (nl1); //&ww[(l*nmean+j)*nfiles];
	      for(i=0;i<=nl;i++) gsl_vector_set (x, i, sqrt(1./nl));
	      gsl_vector_set (x, nl, 0.0);//Good!?
	      T = gsl_multiroot_fsolver_hybrids;
	      s = gsl_multiroot_fsolver_alloc (T, nl1);
	      gsl_multiroot_fsolver_set (s, &df, x);
	      gsl_print_state (iter, s, par);
	      do {
		iter++;
		status = gsl_multiroot_fsolver_iterate (s);
		gsl_print_state (iter, s, par);
		/* check if solver is stuck */
		if (status) {
		  printf ("Minimum not found (ERR=%d)? See /usr/include/gsl/gsl_errno.h for details.",status);
		  if (status==GSL_EBADFUNC) {printf ("Encountered a singular point.\n");}
		  else if (status==GSL_EBADFUNC) {printf ("Noprogress.\n");}
		  else if (status==GSL_ENOPROGJ) {printf (" ERR=%d:Jacobian is not improving the solution.\n",GSL_ENOPROGJ);}
		  break;}//break at Sucess or failure ?
		//		status = gsl_multiroot_test_residual (s->f, 1e-7);
		status = gsl_multiroot_test_residual (s->f, 1./nl);
	      }  while (status == GSL_CONTINUE && iter < itermax);
	      for(i=0;i<nl;i++) ww[(l*nmean+j)*nfiles+i]=square(gsl_vector_get (s->x, i));
	      gsl_multiroot_fsolver_free (s);
	      gsl_vector_free (x);
	      double sumw=0;
	      fprintf(stderr,"w[%d,%d]:",j,l);
	      for(i=0;i<nl;i++) fprintf(stderr,"%+.4f ",ww[(l*nmean+j)*nfiles+i]);
	      for(i=0;i<nl;i++) sumw+=ww[(l*nmean+j)*nfiles+i];
	      fprintf(stderr,"sumw=%f\n",sumw);
	    }
	  }
	}
      }
      else{//for bayesmode!=2
	for(l=0;l<nmean;l++){
	  if(LDmode==2 && l!=m_LD) continue;//only yoko
	  for(j=0;j<nmean;j++){
	    if(LDmode==4 && (l-m_LD)*(j-m_LD)!=0 ) continue;//only tate&yoko
	    for(i=0;i<nfiles2+l-m_LD;i++) ww[(l*nmean+j)*nfiles+i]=1./(nfiles2+l-m_LD);
	  }
	  //	  fprintf(stderr,"w[0-%d,0-%d]=%+.4f\n",nfiles2,nfiles2+l-m_LD,1./(nfiles2+l-m_LD));
	}
      }
      LtrainE=Ltrain=Lvar0=0;
      for(i=0;i<num0;i++){
	//      yp[i]=0;for(j=m_LD;j<nfiles2+m_LD;j++) yp[i]+=(yp01[j*num0+i]*ww[m_LD*nfiles+j-m_LD]);//prediction of training data
	yp[i]=0;
	for(j=0;j<nfiles2;j++){
	  yp[i]+=(yp01[(j+m_LD)*num0+i]*ww0[j]);//prediction of training data
	}
	//	for(j=m_LD;j<nfiles2+m_LD;j++){
	//	  //	  yp[i]+=(yp01[j*num0+i]*ww[(m_LD+m_LD*nmean)*nfiles+j-m_LD]);//prediction of training data
	//	  yp[i]+=(yp01[j*num0+i]*ww0[j-m_LD]);//prediction of training data
	//	}
	yp[i]=y_res(yp[i],&res);
	for(j=m_LD;j<nfiles2+m_LD;j++){
	  Lvar0+= square(yp[i]-yp01[j*num0+i]);
	  LtrainE += square(yp01[(j)*num0+i]-yt[i]);
	}
	//	for(j=m_LD;j<nfiles2+m_LD;j++) Lvar0+= exp(square(yp[i]-yp01[j*num0+i]));//???060530
	Ltrain+= square(yt[i]-yp[i]);
      }
      //      Lvar0 /=(nfiles2*num0); //      Lvar0 = VarY;//?????
      //Lvar0 /=(nfiles2*num0-1+1e-10); //unbiased estimator
      if(nfiles2>=2) Lvar0 /=((num0+ZERO)*(nfiles2)*(nfiles2-1));//estimated variance of ensemble prediction
      else Lvar0 /=(num0+ZERO); //average of unbiased estimator
      Ltrain/=(num0+ZERO);
      LtrainE/=(nfiles*num0+ZERO);
      //      free(yp01);//??why bat result if active?
#define Ltrain4Entropy
#ifdef Ltrain4Entropy
      {
	double pn,E2sum=Ltrain*num0;//Entropy
	H0=0;
	for(i=0;i<num0;i++){
	  pn=square(yp[i]-yt[i])/E2sum;
	  H0 -= pn*log(pn);
	}
	H0/=log(num0+ZERO);
	//	Ltrain=H0;
      }
#endif
      free(yt);//
      free(yp);//
    }//end of "for Lvar0"
#define BAYES2GE 1
#if BAYES2GE == 1
    {//
      int li,l;
      yp00=(double*) malloc(sizeof(double)*num*nfiles);
      //#if LDmode == 2
      //      double *ypbuff=(double*) malloc(sizeof(double)*num*nmean);
      //#elif LDmode == 4
      //      double *ypbuff=(double*) malloc(sizeof(double)*num*nmean*2);
      double *ypbuff=(double*) malloc(sizeof(double)*num*nmean*(nmean));
      //#endif
      y=(double*) malloc(sizeof(double)*num);
      yr=(double*) malloc(sizeof(double)*num);
      n=(int*) malloc(sizeof(int)*num);
      var=(double*) malloc(sizeof(double)*num);
      vv=(double*) malloc(sizeof(double)*num);
      //      cc=(int*) malloc(sizeof(int)*num);
      cc=(int*) malloc(sizeof(int)*num*nfiles);
      kyorix2=(double*)malloc(sizeof(double)*num);
      //      yp=&ypbuff[num*(int)(bayesopt2)];//??what happen when bayesopt2==0??
      yp=&ypbuff[(m_LD+m_LD*nmean)*num];//center of nmean*nmean=(2*m_LD+1)*(2*m_LD+1)
      // yp=&ypbuff[m_LD*num];//
      for(i=0;i<num;i++){yp[i]= n[i]= y[i]= yr[i]=var[i]=vv[i]=0; }
      {///???
	for(j=0;j<nfiles;j++){
	  sprintf(fn_pred,"%spred.dat",argv[j+j0]);
	  fp=fopen(fn_pred,"r");
	  for(i=0;i<num;i++){
	    fgets(buff,buffsize,fp);
	    sscanf(buff,"%lf%d%lf%lf%lf%lf%lf%lf%lf",&_yp,&_n,&_y,&_yr,&_var,&_cc,&_vv,&_e2,&_kyorix2);
	    //	    yp0[i*nfiles+j] =y_res(_yp,&res);
	    if(res.r3<0) yp00[j*num+i]=_yr; //neglect resolution for the meanpred
	    else yp00[j*num+i] =y_res(_yp,&res);
	    //	    yp00[j*num+i] =y_res(_yp,&res);
	    cc[j*num+i] =_cc;//??
	    if(ccmax<_cc) ccmax=_cc;
	    if(j>=m_LD && j<nfiles-m_LD){
	      //	      yp[i]+=_yp*prob[j];//sum for j in [1,...,nfiles-2]
	      n[i]+=_n;
	      y[i]+=_y;//sum for j in [1,...,nfiles-2]
	      yr[i]+=_yr;
	      var[i]+=(_var*_vv);
	      vv[i]+=_vv;
	      //	      cc[i]=_cc;
	      kyorix2[i]=_kyorix2;
	    }
	  }
	  //	  if(j>=m_LD && j<nfiles-m_LD) probtotal+=prob[j];
	  fclose(fp);
	}
	//////////////
	Lvar=Lvarval=Lval=0;
	skew=kurtosis=entropy=mcov=skew2=0,skewa=0;
	m3p=m3m=m2=m3a=m4=m5a=m6=0;nm3p=nm3m=0;
	for(i=0;i<num;i++){//num=1?
	  //	  yp[i]/=probtotal;//mean for j in [1,...,nfiles-2]
	  //	  yp[i]=y_res(yp[i]/probtotal,&res);
	  yp[i]=0;
	  //	  for(j=0;j<nfiles2;j++) {yp[i] +=yp0[i*nfiles+j]*ww[(m_LD)*nfiles+j];}
	  //	  for(j=0;j<nfiles2;j++) {yp[i] +=yp00[(j+m_LD)*num+i]*ww[(m_LD)*nfiles+j];}
	  //	  for(j=m_LD;j<nfiles2+m_LD;j++) yp[i]+=(yp00[j*num+i]*ww[m_LD*nfiles+j]);//prediction of center of test data
	  //	  for(j=0;j<nfiles2;j++) yp[i]+=(yp00[(j+m_LD)*num+i]*ww[(m_LD+n_LD*nmean)*nfiles+j]);//prediction of center of test data
	  for(j=0;j<nfiles2;j++) yp[i]+=(yp00[(j+m_LD)*num+i]*ww0[j]);//prediction of center of test data
	  yp[i]=y_res(yp[i],&res);
	  n[i] = (int)((double)n[i]/nfiles2+0.5);
	  y[i] = y[i]/nfiles2;//mean for j in [1,...,nfiles-2]
	  yr[i]= yr[i]/nfiles2;
	  //	  for(j=m_LD;j<nfiles-m_LD;j++){
	  //	  for(j=m_LD;j<nfiles2+m_LD;j++){
	  double skewj=0,sigma2j=0;
	  double m2j=0,m3j=0,m4j=0,m5j=0,m6j=0;
	  double m3pj=0,m3mj=0;
	  //	  int nm3pj=0,nm3mj=0;
	  double nm3pj=0,nm3mj=0;
	  for(j=0;j<nfiles2;j++){//number of bags
	    double err=yp00[(j+m_LD)*num+i]-yp[i];//inversed 091008
	    double err2_=err*err;
	    m2j+=err2_;//sigma2j+=err2;//m2j
	    skew2+=(err>0)?(err2_):(-err2_);//?
	    m3j+=(err2_*=err);//skewj+=(errn=err2*err);//m3j
	    if(err>0){
	      nm3pj++;
	      m3pj+=err2_;//err^3
	    }
	    else if(err<0){
	      nm3mj++;
	      m3mj+=err2_;//err^3
	    }
	    m4j+=(err2_*=err);//kurtosis+=(errn=errn*err);//m4j
	    m5j+=(err2_*=err);
	    m6j+=(err2_*=err);
	    //Lvar    += square(yp[i]-yp00[(j+m_LD)*num+i]);
	    //Lvar    += square(yp[i]-yp00[(j+m_LD)*num+i]);
	    //Lvar    += exp(square(yp[i]-yp00[(j+m_LD)*num+i]));//???060530
	    Lvarval += square(yp00[(j+m_LD)*num+i]-y[i]);
	    {//mean covariance
	      int l;
	      for(l=0;l<nfiles2;l++){
		if(l!=j) mcov+=(err*(yp00[(l+m_LD)*num+i]-yp[i]));
	      }
	    }
	  }//closing for(j=0;j<nfiles2;j++){
	  Lval += square(yp[i]-y[i]);
	  Lvar+=m2j;//Lvar+=sigma2j;
	  double sj=sqrt(m2j/nfiles2);// double sj=sqrt(sigma2j/nfiles);
	  sj+=ZERO;
	  nm3pj+=ZERO;
	  nm3mj+=ZERO;
	  int mnorm=0;
	  if(mnorm==0){//original
	    m2+=(m2j/nfiles2);
	    m3+=(m3j/nfiles2)/pow(sj,3.);//skew+=skewj;//m3
	    m3a+=fabs(m3j/nfiles2)/pow(sj,3.);//skewa+=fabs(skewj/nfiles2)/pow(sj,3);//m3a
	    m3p+=(m3pj/(nm3pj))/pow(sj,3.);
	    m3m+=(m3mj/(nm3mj))/pow(sj,3.);
	    nm3p+=nm3pj;
	    nm3m+=nm3mj;
	    m4 +=    (m4j/nfiles2)/pow(sj,4.);//
	    m5a+=fabs(m5j/nfiles2)/pow(sj,5.);//?5/2
	    m6 +=    (m6j/nfiles2)/pow(sj,6.);//?3=6/2
	  }
	  else if(mnorm==1){//new for normalize m3,m4,m5,m6
	    m2+=(m2j/nfiles2);
	    m3+=pow(m3j/nfiles2,1./3.)/sj;//skew+=skewj;//m3
	    m3a+=pow(fabs(m3j/nfiles2),1./3.)/sj;//skewa+=fabs(skewj/nfiles2)/pow(sj,3);//m3a
	    m4 +=pow(m4j/nfiles2,1./4.)/sj;//
	    m5a+=pow(fabs(m5j/nfiles2),1./5.)/sj;//?5/2
	    m6 +=pow(m6j/nfiles2,1./6.)/sj;//?3=6/2
	  }
	  else if(mnorm==2){//this does not work for skew and kurtosis
	    m2+=(m2j/nfiles2);
	    m3+=pow(m3j/nfiles2,1./3.);//skew+=skewj;//m3
	    m3a+=pow(fabs(m3j/nfiles2),1./3.);//skewa+=fabs(skewj/nfiles2)/pow(sj,3);//m3a
	    m4 +=pow(m4j/nfiles2,1./4.);//
	    m5a+=pow(fabs(m5j/nfiles2),1./5.);//
	    m6 +=pow(m6j/nfiles2,1./6.);//
	  }
	  //	  if(i<10) fprintf(stderr,"[%d]Lval=%9.3e,yp=%.3e,y=%.3e\n",i,Lval,yp[i],y[i]);
	  //	  if(i>1900) fprintf(stderr,"[%d]Lval=%9.3e,yp=%.3e,y=%.3e\n",i,Lval/num,yp[i],y[i]);
	}//closing for(i=0;i<num;i++){
	int num2=num*nfiles2;
	//Lvar /= (nfiles2*num);
	//??	if(nfiles2>=2) Lvar /=((num)*(nfiles2)*(nfiles2-1.));//estimated variance of ensemble prediction
	//if(nfiles2>=2) Lvar /=num2;//estimated variance of ensemble prediction
	//else Lvar/=(num);
	for(i=0;i<num;i++){
	  for(j=0;j<nfiles2;j++){
	    double err=yp[i]-yp00[(j+m_LD)*num+i];
	    double p=(err*err)/(Lvar+ZERO);
	    p+=ZERO;
	    entropy -= p*log(p);
	  }
	}
	entropy/=log(num2+ZERO);
	//
	Lvar/=num2;//?see above
	mcov/=(num2*(nfiles2-1));
	double sigma=sqrt(Lvar);
	double sigma3=sigma*sigma*sigma;
	m2/=num;
	m3/=num;//skew=(skew/num2)/sigma3;//	skewa/=num;
	m3a/=num;//	skewa/=num;
	m3p/=num;
	m3m/=num;
	m4/=num;//kurtosis=(kurtosis/num2)/sigma3/sigma;//!modified 20130418
	m5a/=num;
	m6=m6/num;
	skew=m3;skewa=m3a;kurtosis=m4;
	Lvar=m2+ZERO;//??
	//	skew2=skew2/num2/Lvar;
	//	kurtosis=(kurtosis/num2)/sigma3/sigma-3;//!!!m4-3
	/////////////
	//Lvar /=((nfiles2*num)*(nfiles2*num-1.));//A-type Uncertainty
	//	Lvar /=((nfiles2*num)*(nfiles2*num-1.));//A-type Uncertainty
	Lvarval /= (num2);
	Lval /= (num);
	//	fprintf(stderr,"#Lvar=%e,Lvarval=%e,Lval=%e,_LD:",Lvar,Lvarval,Lval);
	//	fprintf(stderr,"#Lvar=%e,Lvarval=%e,Lval=%e,num=%d,_LD:",Lvar,Lvarval,Lval,num);
	//fprintf(stderr,"%g %g %g %g %g %g #Ltst,Lvar,skew,krt,ent,log|s*k|#n%db%d\n",Lval,Lvar,skew,kurtosis,entropy,log(fabs(skew*kurtosis)),num,nfiles);
	{
	  FILE *fp=fopen("loss.dat","w");

	  fprintf(fp,"%g %g %g %g %g %d %d %g %g %g %d %g %g %g %g %d %g %d\n#Ltst,Lvar,skw,krt,ent,n,b,cor,skwa,yp0,num,skw2,m5a,m6,m3p,nm3p,m3m,nm3m\n",
		  Lval,Lvar,m3,m4,entropy,num,nfiles,mcov/Lvar,m3a,yp[0],num,skew2,m5a,m6,m3p,nm3p,m3m,nm3m);//4=K,9=sa,13=m5a,14=m6
	  //      1    2    3  4  5       6   7      8         9   10    11  12    13  14 15  16   17  18
	  //	  fprintf(fp,"%g %g %g %g %g %d %d %g %g %g %d %g %g %g\n#Ltst,Lvar,skw,krt,ent,n,b,cor,skwa,yp0,num,skw2\n",Lval,Lvar,skew,kurtosis,entropy,num,nfiles,mcov/Lvar,skewa,yp[0],num,skew2,m5a,m6);//4=K,9=sa,13=m5a,14=m6
	  //	  fprintf(fp,"%g %g %g %g %g %d %d %g %g %g %d\n#Ltst,Lvar,skw,krt,ent,n,b,cor,skew2,yp0,num\n",Lval,Lvar,skew,kurtosis,entropy,num,nfiles,mcov/Lvar,skew2,yp[0],num);
	  //	  fprintf(fp,"%.4g %.4g %.4g %.4g %.4g %d %d\n#Ltst,Lvar,skew,krt,ent,n,b\n",Lval,Lvar,skew,kurtosis,entropy,num,nfiles);
	  fclose(fp);
//	  fprintf(stderr,"###cat loss.dat###\n");
//	  system("cat loss.dat");
//	  fprintf(stderr,"#######\n");
	}
	//	fprintf(stderr,"#Lvar%e,Lvarval%e,Ltst%e,skew%.3g,krt%g,ent%.3g,n%d,b%d,",Lvar,Lvarval,Lval,skew,kurtosis,entropy,num,nfiles);
	int jj;
	for(i=0;i<num;i++){
	  for(j=-m_LD;j<=m_LD;j++){//l1
	    for(l=-m_LD;l<=m_LD;l++){//ld
	      if(j==0 && l==0 ) continue;//already obtained above
	      if(LDmode==2 && l!=0) continue;//only yoko
	      if(LDmode==4 && (l)*(j)!=0 ) continue;//only tate&yoko
	      li=(l*nmean+j)*num+i;
	      yp[li]=0;
	      //	    for(j=0;j<nfiles2;j++) {yp[li] +=yp0[i*nfiles+j]*ww[(m_LD+l)*nfiles+j];}
	      //	    for(j=0;j<nfiles2;j++) {yp[li] +=yp00[(j+m_LD+l)*num+i]*ww[(m_LD+l)*nfiles+j];}
	      for(jj=0;jj<nfiles2+l;jj++) {
		yp[li] +=yp00[(jj+m_LD+j)*num+i]*ww0[(l*nmean+j)*nfiles+jj];
		//		fprintf(stderr,"jj+m_LD+j=%d+%d+%d=%d\n",jj,m_LD,j,jj+m_LD+j);
	      }
	      yp[li]=y_res(yp[li],&res);
	      //	    li=(l+m_LD+1)*num+i;
	      //	    yp[li]=0;
	      //	    for(j=0;j<nfiles2+l;j++) {yp[li] +=yp00[(j+m_LD)*num+i]*ww[((l*nmean+m_LD)*nfiles+j];}
	      //	    yp[li]=y_res(yp[li],&res);
	    }
	  }
	}
	if(1==0){//unnecessary????
	  LD=0;
	  for(j=-m_LD;j<=m_LD;j++){//l1
	    for(l=-m_LD;l<=m_LD;l++){
	      if(j==0 && l==0 ) continue;//already obtained above
	      if(LDmode==2 && l!=0) continue;//only yoko
	      if(LDmode==4 && (l)*(j)!=0 ) continue;//only tate&yoko
	      _LD=0;
	      for(i=0;i<num;i++) {
		_LD += square(yp[i]-yp[(l*nmean+j)*num+i]);
		//		_LD += square(yp[i]-yp[l*num+i]);
		//		_LD += square(yp[i]-yp[(l+m_LD+1)*num+i]);
	      }
	      _LD/=num;
	      LD+=_LD;
	      fprintf(stderr,"_LD%9.3e ",_LD);
	    }
	  }
	  if(LDmode==2) LD/=(m_LD*4.);//
	  else if(LDmode==4) LD/=(m_LD*8.);//
	  else LD/=(nmean*nmean-1);//
	  //	  LD/=(m_LD*8.);//LD/=(m_LD*2.);
	  fprintf(stderr,"LD%e\n",LD);
	}//end of if(1==0)
      }//end of ???
    }//end of for BAYES2GE
#endif
  }//end of bayesmode<100

  //obtain heterosedastic err only vaid for known y[i] ?
#define newvari
#ifdef newvari //20150115
  double *bias=(double*) malloc(sizeof(double)*num);
  double *err=(double*) malloc(sizeof(double)*num);
  int    *nVi=(int*) malloc(sizeof(int)*num);
  {
    int ji,ii;
    for(i=0;i<num;i++) {
      bias[i]=var[i]=nVi[i]=0;
    }
    ccmax+=1;
    double *varcc=(double*)malloc(sizeof(double)*ccmax);
    double *biascc=(double*)malloc(sizeof(double)*ccmax);
    int *nvarcc=(int*)malloc(sizeof(int)*ccmax);
    //    free *varcc;
    for(j=0;j<nfiles;j++){//bags
      int cci;
      for(cci=0;cci<ccmax;cci++) varcc[cci]=biascc[cci]=nvarcc[cci]=0;
      for(i=0;i<num;i++){
	err[i]=yp00[j*num+i]-yp[i];
	varcc[cc[i]]+=square(err[i]);
	biascc[cc[i]]=err[i];
	nvarcc[cc[i]]++;
//	{
//	  fp=fopen("predict+.dat","w");
//	  fprintf(fp,"hello\n");
//	  fclose(fp);
//	}
      }
      for(i=0;i<num;i++){
	var[i]+=varcc[cc[i]];
	bias[i]+=biascc[cc[i]];
	nVi[i]+=nvarcc[cc[i]];
      }
    }
    for(i=0;i<num;i++){
      bias[i]/=nVi[i];
      var[i]/=nVi[i];
    }
    net_printf("##################new variance###\n");
  }

#else //old var[i]
  double *bias=(double*) malloc(sizeof(double)*num);
  double *err=(double*) malloc(sizeof(double)*num);
  int    *nVi=(int*) malloc(sizeof(int)*num);
  int    *nVib=(int*) malloc(sizeof(int)*num);

  {
    int ji,ii;
    for(i=0;i<num;i++) {
      err[i]=yp[i]-y[i];
      bias[i]=var[i]=nViv[i]=nVib[i]=0;
    }
    for(j=0;j<nfiles;j++){//bags
      for(i=0;i<num;i++){
	bias[i]+=err[i];        nVib[i]++;
	var[i] +=square(err[i]);nVi[i]++;
	//	if(1==0){//for check Lval0=9.126e-02
	{//for obtaining the average of the noise
	  ji=j*num+i;
	  for(ii=i+1;ii<num;ii++){
	    if(cc[ji]==cc[j*num+ii]) {//error for the same cell
	      //bias[i]+=err[ii];  nVib[i]++;
	      //bias[ii]+=err[i];  nVib[ii]++;
	      var[i] +=square(err[ii]); nVi[i]++;
	      var[ii]+=square(err[i]);  nVi[ii]++;
	    }
	  }
	}
      }
    }
    for(i=0;i<num;i++){
      bias[i]/=nVib[i];
      var[i]/=nViv[i];
    }
  }
#endif
  //////
  //  fp=fopen(fnpredict,"w");
  fp=fopen("predict+.dat","w");
  for(i=0;i<num;i++){
    if(intmode==1){
      //      fprintf(fp,"%d %d %.7e %.7e %.7e %d %.7e %.7e %.7e #Y^,t,Y,yr,var,c,v,k,b\n",(int)(yp[i]+0.5),n[i],y[i],yr[i],var[i]/vv[i],cc[i],vv[i],kyorix2[i],bias[i]);
      fprintf(fp,"%d %d %.7e %.7e %.7e %d %.7e %.7e %.7e %d #Y^,t,Y,yr,var,c,v,k,b,nvar\n",(int)(yp[i]+0.5),n[i],y[i],yr[i],var[i],cc[i],vv[i],kyorix2[i],bias[i],nVi[i]);
    }
    else{
      fprintf(fp,"%.7e %d %.7e %.7e %.7e %d %.7e %.7e %.7e %d #Y^,t,Y,yr,var,c,v,k,b,nvar\n",yp[i],n[i],y[i],yr[i],var[i],cc[i],vv[i],kyorix2[i],bias[i],nVi[i]);
//      {//error check
//	char buff[256];
//	sprintf(buff,"%.7e %d %.7e %.7e %.7e %d %.7e %.7e %.7e %d #Y^,t,Y,yr,var,c,v,k,b,nvar\n",yp[i],n[i],y[i],yr[i],var[i],cc[i],vv[i],kyorix2[i],bias[i],nVi[i]);
//	fprintf(stdout,"%s",buff);
//	fprintf(fp,"%s",buff);
//      }
      //      fprintf(fp,"%.7e %d %.7e %.7e %.7e %d %.7e %.7e %.7e #Y^,t,Y,yr,var,c,v,k,b\n",yp[i],n[i],y[i],yr[i],var[i]/vv[i],cc[i],vv[i],kyorix2[i],bias[i]);
    }
  }
  fclose(fp);
  //  fp=fopen(fnloss,"w");
  fp=fopen("loss+.dat","w");
  fprintf(fp,"%13.7e %13.7e %13.7e %13.7e %13.7e %d %d %.7e %.7e %.7e #LD Lval Lvar Lvarval Lvar0 nfiles num Ltrain LtrainE H0\n",LD,Lval,Lvar,Lvarval,Lvar0,nfiles,num,Ltrain,LtrainE,H0);
  fclose(fp);
  //  fprintf(stderr,"Created '%s'=%s",fnpredict,fnbody(argv[1],pr));
  //  for(j=2;j<=nfiles;j++) fprintf(stderr,"+%s" ,fnbody(argv[j],pr));
  //  fprintf(stderr,"\n");
  net_printf(strx3("Created '%s' and '%s'. #LDmode=%d.\n","predict+.dat","loss+.dat",LDmode));
  return(1);
}

ext: c


text:
#include <math.h>
#include <stdlib.h>

#include <fastpm/libfastpm.h>
#include <fastpm/logging.h>
#include <fastpm/constrainedgaussian.h>

#include "pmpfft.h"
#include <gsl/gsl_linalg.h>

double fastpm_2pcf_eval(FastPM2PCF* self, double r)
{
    double rMax = self->size * self->step_size;
    if(r > rMax)
        return 0.0;
    if(r == rMax)
        return self->xi[self->size];
    int i = (int)floor(r / self->step_size);
    double prev = self->xi[i], next = self->xi[i + 1];
    double deltaR = r - i * self->step_size;
    return prev + (next - prev) * deltaR / self->step_size;
}

/*
void
fastpm_generate_covariance_matrix(PM *pm, fastpm_fkfunc pkfunc, void * data, FastPMFloat *cov_x)
{
}
*/
void
fastpm_2pcf_from_powerspectrum(FastPM2PCF *self, fastpm_fkfunc pkfunc, void * data, double r_max, int steps)
{
    self->size = steps;
    self->step_size = r_max / steps;
    self->xi = malloc((steps + 1) * sizeof(double));

    double logKMin = -10, logKMax = 5;
    double logKSteps = 10000;
    double logKStepSize = (logKMax - logKMin) / logKSteps;
    double pi = 3.141593;
    int i, j;
    for(i = 0; i <= steps; ++i)
    {
        double r = i * self->step_size;
        double res = 0;
        double prev = 0;
        for(j = 1; j <= logKSteps; ++j)
        {
            double logK = logKMin + j * logKStepSize;
            double k = exp(logK);
            double kr = k * r;
            double func = 1;
            if(kr > 0)
                func = sin(kr) / kr;
            func *= pkfunc(k, data) * k * k * k;
            res += (prev + func) / 2;
            prev = func;
        }
        res *= logKStepSize / (2 * pi * pi);
        self->xi[i] = res;
    }
}

static void
_solve(int size, double * Cij, double * dfi, double * x)
{
    gsl_matrix_view m = gsl_matrix_view_array(Cij, size, size);
    gsl_vector_view b = gsl_vector_view_array(dfi, size);
    gsl_vector_view xx = gsl_vector_view_array(x, size);
    gsl_permutation *p = gsl_permutation_alloc(size);
    int s;
    gsl_linalg_LU_decomp(&m.matrix, p, &s);
    gsl_linalg_LU_solve(&m.matrix, p, &b.vector, &xx.vector);
}

static void
_readout(FastPMConstraint * constraints, int size, PM * pm, FastPMFloat * delta_x, double * dfi)
{
    int i;

    for(i = 0; i < size; ++i)
    {
        int ii[3];
        int d;
        int inBox = 1;
        int index = 0;
        for(d = 0; d < 3; ++d)
        {
            ii[d] = constraints[i].x[d] * pm->InvCellSize[d] - pm->IRegion.start[d];
            if(ii[d] < 0 || ii[d] > pm->IRegion.size[d])
                inBox = 0;
            index += ii[d] * pm->IRegion.strides[d];
        }

        if(inBox) {
            dfi[i] = delta_x[index];
        } else {
            dfi[i] = 0;
        }
    }
    MPI_Allreduce(MPI_IN_PLACE, dfi, size, MPI_DOUBLE, MPI_SUM, pm_comm(pm));
}
static double
_sigma(PM * pm, FastPMFloat * delta_x)
{
    double d2 = 0.0;
    PMXIter xiter;
    for(pm_xiter_init(pm, &xiter);
       !pm_xiter_stop(&xiter);
        pm_xiter_next(&xiter))
    {
        double od = delta_x[xiter.ind] - 1;
        d2 += od * od;
    }
    MPI_Allreduce(MPI_IN_PLACE, &d2, 1, MPI_DOUBLE, MPI_SUM, pm_comm(pm));
    /* unbiased estimator of the variance. the mean is 1. */
    d2 /= (pm_norm(pm) - 1);
    return sqrt(d2);
}

void
fastpm_cg_apply_constraints(FastPMConstrainedGaussian *cg, PM * pm, FastPM2PCF *xi, FastPMFloat * delta_k)
{
    FastPMConstraint * constraints = cg->constraints;

    int size;
    for(size = 0; constraints[size].x[0] >= 0; size ++)
        continue;

    int i;
    fastpm_info("Constrained Gaussian with %d constraints\n", size);
    double dfi[size];
    double e[size];
    double Cij[size * size];
    double sigma = 0;
    FastPMFloat * delta_x = pm_alloc(pm);

    pm_assign(pm, delta_k, delta_x);
    pm_c2r(pm, delta_x);

    sigma = _sigma(pm, delta_x);
    fastpm_info("Measured sigma on the grid = %g\n", sigma);

    _readout(constraints, size, pm, delta_x, dfi);

    for(i = 0; i < size; ++i)
    {
        dfi[i] = (1 + constraints[i].c * sigma) - dfi[i];

        int j;
        for(j = i; j < size; ++j)
        {
            int d;
            double r = 0;
            for(d = 0; d < 3; ++d) {
                double dx = constraints[i].x[d] - constraints[j].x[d];

                if(dx > 0.5*pm->BoxSize[d]){
                    dx -= pm->BoxSize[d];
                }
                else if(dx < -0.5*pm->BoxSize[d]){
                    dx += pm->BoxSize[d];
                }

                r += dx * dx;
            }
            r = sqrt(r);
            double v = fastpm_2pcf_eval(xi, r);
            Cij[i * size + j] = v;
            Cij[j * size + i] = v;
        }
    }

    _solve(size, Cij, dfi, e);

    PMXIter xiter;
    for(pm_xiter_init(pm, &xiter);
       !pm_xiter_stop(&xiter);
        pm_xiter_next(&xiter))
    {
        double v = 0;
        for(i = 0; i < size; ++i)
        {
            int d;

            double r = 0;
            for(d = 0; d < 3; d ++) {
                double dx = xiter.iabs[d] * pm->CellSize[d] - constraints[i].x[d];

                if(dx > 0.5*pm->BoxSize[d]){
                    dx -= pm->BoxSize[d];
                }
                else if(dx < -0.5*pm->BoxSize[d]){
                    dx += pm->BoxSize[d];
                }

                r += dx * dx;
            }
            r = sqrt(r);

            v += e[i] * fastpm_2pcf_eval(xi, r);
        }
        delta_x[xiter.ind] += v;
    }

    _readout(constraints, size, pm, delta_x, dfi);
    for(i = 0; i < size; i ++) {
        fastpm_info("After constraints, Realization x[] = %g %g %g overdensity = %g, peak-sigma= %g\n",
                constraints[i].x[0],
                constraints[i].x[1],
                constraints[i].x[2],
                (dfi[i] - 1.0), (dfi[i] - 1.0) / sigma);
    }
    pm_r2c(pm, delta_x, delta_k);
    pm_free(pm, delta_x);
}

max_stars_repo_name: sbird/FastPMRunner

max_issues_repo_name: sbird/FastPMRunner

max_forks_repo_name: sbird/FastPMRunner


text:
#include <stdlib.h>
#include <math.h>

#include <fftw3.h>

#include <gsl/gsl_sf_result.h>
#include <gsl/gsl_sf_gamma.h>

#include "fftlog.h"

#ifndef M_PI
    #define M_PI 3.14159265358979323846
#endif

/* This code is FFTLog, which is described in arXiv:astro-ph/9905191 */

static double_complex lngamma_fftlog(double_complex z)
{
  gsl_sf_result lnr, phi;
  gsl_sf_lngamma_complex_e(creal(z), cimag(z), &lnr, &phi);
#ifdef _COMPLEX_STRUCT
  double_complex res = {lnr.val, phi.val};
  return res;
#else
  return lnr.val + I*phi.val;
#endif
}

static double_complex gamma_fftlog(double_complex z)
{
  return cexp(lngamma_fftlog(z));
}

static double_complex polar (double r, double phi)
{
#ifdef _COMPLEX_STRUCT
  double_complex res = {r*cos(phi) , r*sin(phi)};
  return res;
#else
  return (r*cos(phi) +I*(r*sin(phi)));
#endif
}

static void lngamma_4(double x, double y, double* lnr, double* arg)
{
#ifdef _COMPLEX_STRUCT
  double_complex z = {x,y};
  double_complex w = lngamma_fftlog(z);
#else
  double complex w = lngamma_fftlog(x+y*I);
#endif
  if(lnr) *lnr = creal(w);
  if(arg) *arg = cimag(w);
}

static double goodkr(int N, double mu, double q, double L, double kr)
{
  double xp = (mu+1+q)/2;
  double xm = (mu+1-q)/2;
  double y = M_PI*N/(2*L);
  double lnr, argm, argp;
  lngamma_4(xp, y, &lnr, &argp);
  lngamma_4(xm, y, &lnr, &argm);
  double arg = log(2/kr) * N/L + (argp + argm)/M_PI;
  double iarg = round(arg);
  if(arg != iarg)
    kr *= exp((arg - iarg)*L/N);
  return kr;
}

void compute_u_coefficients(int N, double mu, double q, double L, double kcrc, double_complex u[])
{
  double y = M_PI/L;
  double k0r0 = kcrc * exp(-L);
  double t = -2*y*log(k0r0/2);

  if(q == 0) {
    double x = (mu+1)/2;
    double lnr, phi;
    for(int m = 0; m <= N/2; m++) {
      lngamma_4(x, m*y, &lnr, &phi);
      u[m] = polar(1.0,m*t + 2*phi);
    }
  }
  else {
    double xp = (mu+1+q)/2;
    double xm = (mu+1-q)/2;
    double lnrp, phip, lnrm, phim;
    for(int m = 0; m <= N/2; m++) {
      lngamma_4(xp, m*y, &lnrp, &phip);
      lngamma_4(xm, m*y, &lnrm, &phim);
      u[m] = polar(exp(q*log(2) + lnrp - lnrm), m*t + phip - phim);
    }
  }

  for(int m = N/2+1; m < N; m++)
    u[m] = conj(u[N-m]);
  if((N % 2) == 0) {
#ifdef _COMPLEX_STRUCT
    double_complex z = {creal(u[N/2]) , 0.0};
    u[N/2] =  z;
#else
    u[N/2] = (creal(u[N/2]) + I*0.0);
#endif
   }
}

void fht(int N, const double r[], const double_complex a[], double k[], double_complex b[], double mu,
         double q, double kcrc, int noring, double_complex* u)
{
  double L = log(r[N-1]/r[0]) * N/(N-1.);
  double_complex* ulocal = NULL;
  if(u == NULL) {
    if(noring)
      kcrc = goodkr(N, mu, q, L, kcrc);
    ulocal = malloc (sizeof(double_complex)*N);
    compute_u_coefficients(N, mu, q, L, kcrc, ulocal);
    u = ulocal;
  }

  /* Compute the convolution b = a*u using FFTs */
  fftw_plan forward_plan = fftw_plan_dft_1d(N, (fftw_complex*) a, (fftw_complex*) b,  -1, FFTW_ESTIMATE);
  fftw_plan reverse_plan = fftw_plan_dft_1d(N, (fftw_complex*) b, (fftw_complex*) b, +1, FFTW_ESTIMATE);
  fftw_execute(forward_plan);
  for(int m = 0; m < N; m++)
#ifdef _COMPLEX_STRUCT
  { double_complex z = _Cmulcc(b[m],u[m]);
    z = _Cmulcr(z, 1/(double)(N));
    b[m] = z;
  }
#else
    b[m] *= u[m] / (double)(N);       // divide by N since FFTW doesn't normalize the inverse FFT
#endif
  fftw_execute(reverse_plan);
  fftw_destroy_plan(forward_plan);
  fftw_destroy_plan(reverse_plan);

  /* Reverse b array */
  double_complex tmp;
  for(int n = 0; n < N/2; n++) {
    tmp = b[n];
    b[n] = b[N-n-1];
    b[N-n-1] = tmp;
  }

  /* Compute k's corresponding to input r's */
  double k0r0 = kcrc * exp(-L);
  k[0] = k0r0/r[0];
  for(int n = 1; n < N; n++)
    k[n] = k[0] * exp(n*L/N);

  free(ulocal);
}

void fftlog_ComputeXi2D(double bessel_order,int N,const double l[],const double cl[],
			double th[], double xi[])
{
  double_complex* a = malloc(sizeof(double_complex)*N);
  double_complex* b = malloc(sizeof(double_complex)*N);

  for(int i=0;i<N;i++){
#ifdef _COMPLEX_STRUCT
    double_complex z =  {l[i]*cl[i],0.0};
    a[i]=z;
#else
	a[i]=l[i]*cl[i];
#endif
    }
  fht(N,l,a,th,b,bessel_order,0,1,1,NULL);
  for(int i=0;i<N;i++)
    xi[i]=creal(b[i])/(2*M_PI*th[i]);

  free(a);
  free(b);
}

void fftlog_ComputeXiLM(double l, double m, int N, const double k[], const double pk[],
			double r[], double xi[])
{
  double_complex* a = malloc(sizeof(double_complex)*N);
  double_complex* b = malloc(sizeof(double_complex)*N);

  for(int i = 0; i < N; i++){
#ifdef _COMPLEX_STRUCT
    double_complex z= { pow(k[i], m - 0.5) * pk[i],0};
    a[i] = z;
#else
	a[i] = pow(k[i], m - 0.5) * pk[i];
#endif
  }
  fht(N, k, a, r, b, l + 0.5, 0, 1, 1, NULL);
  for(int i = 0; i < N; i++)
    xi[i] = pow(2*M_PI*r[i], -(m-0.5)) * creal(b[i]);

  free(a);
  free(b);
}

void pk2xi(int N, const double k[], const double pk[], double r[], double xi[])
{
  fftlog_ComputeXiLM(0, 2, N, k, pk, r, xi);
}

void xi2pk(int N, const double r[], const double xi[], double k[], double pk[])
{
  static const double TwoPiCubed = 8*M_PI*M_PI*M_PI;
  fftlog_ComputeXiLM(0, 2, N, r, xi, k, pk);
  for(int j = 0; j < N; j++)
    pk[j] *= TwoPiCubed;
}

max_stars_repo_name: cmbant/CCL

max_issues_repo_name: cmbant/CCL

max_forks_repo_name: cmbant/CCL


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

#include <gsl/gsl_errno.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_complex.h>

#include "constants.h"
#include "struct.h"
#include "fresnel.h"

#define acc 1e-8

/* Used for previous version of branching between cases 1a, 1b and 2 in ComputeInt */
/* static const double betathreshold = 1./10; */
/* static const double invbetathreshold = 10.; */
/* static const double lambdathreshold = 0.01; */

/* Used for previous version of branching between cases 1a, 1b and 2 in ComputeInt */
/* static const double factorbranchingcondition = 0.2154434690031884; */
/* static const double case1athreshold = 0.1; */
/* Used for case 0 - hoping to cover previous instability */
/* static const double case0threshold = 0.05; */

/* Used for branching between cases 1a, 1b, 2, 3 and 4 in ComputeInt */
/* FIXME: another round of refinement would be welcome, case 3 has some cases with 1e-4 error */
static const double p2threshold = 0.02;
static const double p1threshold1 = 0.1;
static const double p1threshold2 = 100.;
static const double p2p1slopethreshold = 0.2;

static const double coeffdeltaC[13] = {-0.1,-0.0462962962963,-0.0230769230769,-0.0136554621849,-0.00899470899471,-0.00636363636364,-0.00473664266768,-0.00366161616162,-0.00291467938527,-0.00237483953787,-0.0019721019721,-0.00166370896185,-0.00142235123367};
static const double coeffdeltaS[13] = {-0.0714285714286,-0.0318181818182,-0.0174603174603,-0.0109649122807,-0.00750988142292,-0.00546058879392,-0.00414746543779,-0.00325630252101,-0.00262408157145,-0.00215946843854,-0.00180809015222,-0.00153594771242,-0.0013209013209};

static const double fcoeff[12] = {0.3989422275318915, 1.839999368141345e-7, -0.2992429124371606, 0.0029447692036381, 2.48395059327971, 3.897790465813753, -140.2709114924737, 952.371347498219, -3660.703989590287, 8649.34061360465, -11740.93693284928, 7032.8915074629};
static const double gcoeff[12] = {0, 0.1994718004197629, -0.000125952146250442, -0.7386823670326993, -0.353224701618227, 19.45013452321292, -97.8846433217499, 221.1470704754588, -32.29876222016862, -1102.56036216157, 2554.184335427029, -1962.422566478277};

static double invn[12] = {1., 0.5, 0.3333333333333333, 0.25, 0.2, 0.1666666666666667, 0.1428571428571428, 0.125, 0.1111111111111111, 0.1, 0.0909090909090909, 0.0833333333333333};

/* Use for Chebyshev, case 3 */
static double ChebyshevNodes[6] = {0.965925826289068,0.7071067811865475,0.2588190451025207,-0.2588190451025207,-0.7071067811865475,-0.965925826289068};
static double ChebyshevWeights[6][6] = {{0.3333333333333333,0.3333333333333333,0.3333333333333333,0.3333333333333333,0.3333333333333333,0.3333333333333333},{0.3219752754296894,0.2357022603955158,0.0862730150341736,-0.0862730150341736,-0.2357022603955158,-0.3219752754296894},{0.2886751345948129,0.,-0.2886751345948129,-0.2886751345948129,0.,0.2886751345948129},{0.2357022603955158,-0.2357022603955158,-0.2357022603955158,0.2357022603955158,0.2357022603955158,-0.2357022603955158},{0.1666666666666667,-0.3333333333333333,0.1666666666666667,0.1666666666666667,-0.3333333333333333,0.1666666666666667},{0.0862730150341736,-0.2357022603955158,0.3219752754296894,-0.3219752754296894,0.2357022603955158,-0.0862730150341736}};

/* Function computing cexp(I*x)-1 without loss of accuracy*/
static double complex cexpm1i(double x){
  double rr;//rr=cos(x)-1=2*sin(x/2)**2
  double ri=sin(x);
  double sinhalf=sin(0.5*x);
  rr=-2.0*sinhalf*sinhalf;
  return rr+I*ri;
};

/* Function computing the Fresnel integral for any real (positive) number, with 1e-8 accuracy */
/* Reference used: Computation of Fresnel Integrals. II by K. Mielenz, [J. Res. Natl. Inst. Stand. Technol. 105, 589 (2000)] */
static double complex FresnelE(const double x) {
  double complex res;

  /* Check argument */
  if(x<0) return -FresnelE(-x);

  /* If x<2, Taylor expansion */
  else if(x<2) {
    double C = 0;
    double S = 0;
    double x3 = x*x*x;
    double x4 = x3*x;
    int n = 0;
    double deltaC = x;
    double deltaS = x3/3.;
    while(true) {
      C += deltaC;
      S += deltaS;
      //deltaC = -(4*n+1.)/(4*n+5.)/(2*n+1.)/(2*n+2.) * x4 * deltaC;
      //deltaS = -(4*n+3.)/(4*n+7.)/(2*n+2.)/(2*n+3.) * x4 * deltaS;
      deltaC = coeffdeltaC[n] * x4 * deltaC;
      deltaS = coeffdeltaS[n] * x4 * deltaS;
      n++;
      if( (fabs(deltaC)<acc && fabs(deltaS)<acc) || n>12) break;
    }
    res = C + I*S;
  }

  /* If x>2, Mielenz-Boersma formula using tabulated coeffs */
  else {
    double xinv = 1./x;
    double xinv2 = xinv*xinv;
    double y[12] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
    y[0] = xinv;
    for(int i=1; i<12; i++) y[i] = y[i-1]*xinv2;
    double f = 0;
    double g = 0;
    for(int i=0; i<12; i++) {
      f += fcoeff[i] * y[i];
      g += gcoeff[i] * y[i];
    }
    res = sqrt(PI/2) * ( (1+I)/2 - (g+I*f) * cexp(I*x*x));
  }

  return res;
}

/* In fact equivalent to Case 2 */
//static double complex ComputeIntCase0(
//   gsl_vector* coeffsAreal,         /* */
//   gsl_vector* coeffsAimag,         /* */
//   double p0,                       /* */
//   double p1,                       /* */
//   double p2,                       /* */
//   double eps)                      /* */
// {
//   /* Get poly coeffs and rescale by epsilon */
//   double complex coeffs[4] = {0.,0.,0.,0.};
//   double epspow[4] = {0.,0.,0.,0.};
//   epspow[0] = 1.; epspow[1] = eps; epspow[2] = eps*eps; epspow[3] = eps * epspow[2];
//   for(int i=0; i<=3; i++) {
//     coeffs[i] = epspow[i] * (gsl_vector_get(coeffsAreal, i+1) + I*gsl_vector_get(coeffsAimag, i+1));
//   }
//   double p0r = p0;
//   double p1r = p1*epspow[1];
//   double p2r = p2*epspow[2];
//
//   /* In case p2<0, conjugation */
//   int conjug = 0;
//   if(p2r<0) {
//     conjug=1;
//     p2r = -p2r;
//     p1r = -p1r;
//     p0r = -p0r;
//     for(int i=0; i<=3; i++) coeffs[i] = conj(coeffs[i]) ;
//   }
//
//   //
//   //printf("case0\n");
//
//   /* Prefactor and bounds after change of variable */
//   double sqrtp2r = sqrt(p2r);
//   double complex factor = 1./sqrtp2r*cexp(I*p0r);
//   double a = p1r/(2*sqrtp2r);
//   double b = p1r/(2*sqrtp2r) + sqrtp2r;
//   /* Coeffs in polynomial changed by shift */
//   double alpha = 1./sqrtp2r;
//   double beta = -p1r/(2*p2r);
//   double complex a0 = coeffs[0]; double complex a1 = coeffs[1]; double complex a2 = coeffs[2]; double complex a3 = coeffs[3];
//   double alpha2 = alpha*alpha; double beta2 = beta*beta;
//   double complex a0p = a0 + a1*beta + a2*beta2 + a3*beta2*beta;
//   double complex a1p = a1*alpha + 2*a2*alpha*beta + 3*a3*alpha*beta2;
//   double complex a2p = a2*alpha2 + 3*a3*alpha2*beta;
//   double complex a3p = a3*alpha2*alpha;
//
//   /* Result */
//   double complex res = 0.;
//   res = factor*((a0p + I/2*a2p) * (FresnelE(b) - FresnelE(a)) + (1./2*(a3p - I*a1p) - I/2*a2p*b - I/2*a3p*b*b) * cexp(I*b*b) - (1./2*(a3p - I*a1p) - I/2*a2p*a - I/2*a3p*a*a) * cexp(I*a*a));
//   //
//   //printf("%g, %g, %g \n", creal(factor), creal((a0p + I/2*a2p) * (FresnelE(b) - FresnelE(a))), creal(+ (1./2*(a3p - I*a1p) - I/2*a2p*b - I/2*a3p*b*b) * cexp(I*b*b) - (1./2*(a3p - I*a1p) - I/2*a2p*a - I/2*a3p*a*a) * cexp(I*a*a)));
//   //exit(0);
//
//   /* Remembering conjugation */
//   if(conjug) res = conj(res);
//   return res;
// }

double complex ComputeIntCase1a(
  const double complex* coeffsA,         /* */
  const double p1,                       /* */
  const double p2,                       /* */
  const double scale)                    /* */
{
  /* Setting up the coefficients and the required accuracy */
  double complex coeffs[12] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
  for(int i=0; i<=3; i++) coeffs[i] = coeffsA[i];
  int maxj0 = 3;
  double acctol = 1.e-5/scale;

  /* Polynomial given by the expansion of cexp(I*p1*x) */
  double complex poly1[8] = {0.,0.,0.,0.,0.,0.,0.,0.};
  double p1abs = fabs(p1);
  double complex term1 = 1.; double term1absmax = 1.;
  int maxj1 = 0;
  for(int i=0; i<8; i++) {
    poly1[i] = term1;
    term1absmax = term1absmax * p1abs * invn[i]; /* Note: invn[i] = 1/(i+1) */
    if(term1absmax<acctol) break;
    term1 = term1 * I*p1 * invn[i];
    maxj1++;
  }

  /* Polynomial given by the expansion of cexp(I*p2*x2) */
  double complex poly2[11] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
  double p2abs = fabs(p2);
  double complex term2 = 1.; double term2absmax = 1.;
  int maxj2 = 0;
  for(int i=0; i<4; i++) {
    poly2[2*i] = term2;
    term2absmax = term2absmax * p2abs * invn[i]; /* Note: invn[i] = 1/(i+1) */
    if(term2absmax<acctol) break;
    term2 = term2 * I*p2 * invn[i];
    maxj2 += 2;
  }

  /* Multiplying in place the polynomials - we cut according to the fact that we allow for only 12 coeffs up to x^11 (assuming the rest would be negligible anyway), and we know the constant in poly1 and poly2 is always 1. */
  for(int i=min(11, maxj0+maxj1); i>=0; i--) {
    for(int j=max(1, i-maxj0); j<=min(i, maxj1); j++) {
      coeffs[i] += coeffs[i-j] * poly1[j];
    }
  }
  maxj0 = min(11, maxj0 + maxj1);
  for(int i=min(11, maxj0+maxj2); i>=0; i--) {
    for(int j=max(1, i-maxj0); j<=min(i, maxj2); j++) {
      coeffs[i] += coeffs[i-j] * poly2[j];
    }
  }
  maxj0 = min(11, maxj0 + maxj2);

  /* Computing the integral itself */
  double complex res = 0.;
  for(int i=0; i<=maxj0; i++) {
    res += coeffs[i] * invn[i]; /* Note: invn[i] = 1/(i+1) */
  }
  return res;
}

double complex ComputeIntCase1b(
  const double complex* coeffsA,         /* */
  const double p1,                       /* */
  const double p2,                       /* */
  const double scale)                    /* */
{
  /* Setting up the coefficients and the required accuracy */
  double complex coeffs[12] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
  for(int i=0; i<=3; i++) coeffs[i] = coeffsA[i];
  int maxj0 = 3;
  double acctol = 1.e-5/scale;

  /* Polynomial given by the expansion of cexp(I*p2*x2) */
  double complex poly2[11] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
  double p2abs = fabs(p2);
  double complex term2 = 1.; double term2absmax = 1.;
  int maxj2 = 0;
  for(int i=0; i<4; i++) {
    poly2[2*i] = term2;
    term2absmax = term2absmax * p2abs * invn[i]; /* Note: invn[i] = 1/(i+1) */
    if(term2absmax<acctol) break;
    term2 = term2 * I*p2 * invn[i];
    maxj2 += 2;
  }

  /* Multiplying in place the polynomials - we cut according to the fact that we allow for only 12 coeffs up to x^11 (assuming the rest would be negligible anyway), and we know the constant in poly2 is always 1. */
  for(int i=min(11, maxj0+maxj2); i>=0; i--) {
    for(int j=max(1, i-maxj0); j<=min(i, maxj2); j++) {
      coeffs[i] += coeffs[i-j] * poly2[j];
    }
  }
  maxj0 = min(11, maxj0 + maxj2);

  /* Go down the recursion relation */
  double invp1 = 1./p1;
  double complex coeffE1 = 0.; double complex coeffE1minus1 = 0;
  for(int i=maxj0; i>=1; i--) {
    coeffE1 += -I*invp1 * coeffs[i];
    coeffs[i-1] += I*invp1 * i * coeffs[i];
  }
  coeffE1minus1 = -I*invp1 * coeffs[0];
  /* Result */
  double complex cexpm1ip1 = cexpm1i(p1);
  double complex res = (coeffE1*(cexpm1ip1+1.) + coeffE1minus1*cexpm1ip1);
  return res;
}

double complex ComputeIntCase2(
  const double complex* coeffsA,         /* */
  const double p1,                       /* */
  const double p2)                       /* */
{
  /* Setting up the coefficients */
  double complex coeffs[4] = {0.,0.,0.,0.};
  for(int i=0; i<=3; i++) coeffs[i] = coeffsA[i];
  double complex coeffE2 = 0; double complex coefffresnel = 0; double complex constant = 0;
  /* Conjugate if required */
  int conjug = 0; double p1b = p1; double p2b = p2;
  if(p2<0) { conjug = 1; p1b = -p1; p2b = -p2; for(int i=0; i<=3; i++) coeffs[i] = conj(coeffs[i]); }

  /* n=3 and n=2 */
  double inv2p2 = 1. / (2.*p2b);
  for(int i=3; i>=2; i--) {
    coeffE2 += -I*inv2p2 * coeffs[i];
    coeffs[i-2] += I*inv2p2 * (i-1) * coeffs[i];
    coeffs[i-1] += -p1b * inv2p2 * coeffs[i];
  }
  /* n=1 */
  coeffE2 += -I*inv2p2 * coeffs[1];
  constant = I*inv2p2 * coeffs[1];
  coefffresnel += -p1b * inv2p2 * coeffs[1];
  /* n=0 */
  coefffresnel += coeffs[0];

  /* Fresnel integral and complex exponential */
  double complex E2 = cexp(I * (p1b + p2b));
  double sqrtp2 = sqrt(p2b);
  double invsqrtp2 = 1./sqrtp2;
  double halfratio = 0.5 * p1b * invsqrtp2;
  double complex fresnel = cexp(-I*halfratio*halfratio) * invsqrtp2 * (FresnelE(halfratio + sqrtp2) - FresnelE(halfratio));

  /* Result */
  double complex res = (constant + coeffE2*E2 + coefffresnel*fresnel);
  if(conjug) res = conj(res); /* Remember possible conjugation */
  return res;
}

double complex ComputeIntCase3(
  const double complex* coeffsA,         /* */
  const double p1,                       /* */
  const double p2)                       /* */
{
  /* Prepare change of variables */
  double r = p2/p1;

  /* Compute integrand F on nodes */
  double Fvalues[6] = {0.,0.,0.,0.,0.,0.};
  double xk; double y; double sqrtterm; double factor; double x;
  for(int i=0; i<6; i++) {
    xk = ChebyshevNodes[i];
    y = (1.+xk)/2.;
    sqrtterm = sqrt(1. + 4.*r*(1.+r)*y);
    factor = (1.+r)/sqrtterm;
    x = (2.*(1.+r)*y)/(1. + sqrtterm);
    Fvalues[i] = coeffsA[3];
    for(int j=2; j>=0; j--) Fvalues[i] = coeffsA[j] + x*Fvalues[i];
    Fvalues[i] = factor * Fvalues[i];
  }

  /* Compute Chebyshev coefficients ci */
  double complex c[6] = {0.,0.,0.,0.,0.,0.};
  for(int i=0; i<6; i++) {
    for(int j=0; j<6; j++) c[i] += ChebyshevWeights[i][j] * Fvalues[j];
  }

  /* Coefficients of polynomial in 1/(I*p), p=p1*(1+r) */
  double p = p1*(1.+r);
  double complex coeffszcos[6] = {0.,0.,0.,0.,0.,0.};
  double complex coeffszsin[6] = {0.,0.,0.,0.,0.,0.};
  coeffszcos[0] = c[1] + c[3] + c[5];
  coeffszcos[1] = -8.*c[2] - 32.*c[4];
  coeffszcos[2] = 96.*c[3] + 800.*c[5];
  coeffszcos[3] = -1536.*c[4];
  coeffszcos[4] = 30720.*c[5];
  coeffszcos[5] = 0.;
  coeffszsin[0] = 0.5*c[0] + c[2] + c[4];
  coeffszsin[1] = -2.*c[1] - 18.*c[3] - 50.*c[5];
  coeffszsin[2] = 16.*c[2] + 320.*c[4];
  coeffszsin[3] = -192.*c[3] -6720.*c[5];
  coeffszsin[4] = 3072.*c[4];
  coeffszsin[5] = -61440.*c[5];

  /* Coefficients of cos(p/2) and sin(p/2) */
  double complex z = 1./(I*p);
  double complex coeffcos = coeffszcos[5];
  for(int i=4; i>=0; i--) coeffcos = coeffszcos[i] + z*coeffcos;
  double complex coeffsin = coeffszsin[5];
  for(int i=4; i>=0; i--) coeffsin = coeffszsin[i] + z*coeffsin;
  coeffcos = z * coeffcos;
  coeffsin = z * coeffsin;

  /* Result */
  double pov2 = p/2.;
  return cexp(I*pov2) * (2*cos(pov2)*coeffcos + 2*I*sin(pov2)*coeffsin);
}

double complex ComputeIntCase4(
  const double complex* coeffsA,         /* */
  const double p1,                       /* */
  const double p2)                       /* */
{
  double complex a0 = coeffsA[0];
  double complex a1 = coeffsA[1];
  double complex a2 = coeffsA[2];
  double complex a3 = coeffsA[3];
  double r = 2.*p2/p1;
  double inv1plusr = 1./(1.+r); double inv1plusr2 = inv1plusr*inv1plusr;
  double inv1plusr4 = inv1plusr2*inv1plusr2;
  double invp1 = 1./p1; double invp12 = invp1*invp1; double invp13 = invp12*invp1;

  double complex term1 = cexp(I*(p1 + p2))*(-I*invp1*(a0 + a1 + a2 + a3)*inv1plusr + invp12*(a1 + 2.*a2 + 3.*a3)*inv1plusr2 + I*invp13*(2.*a2 + 6.*a3 + r*(3.*a3 - a1))*inv1plusr4);
  double complex term2 = - (-I*invp1*a0 + invp12*a1 + I*invp13*(2.*a2 - a1*r));
  return term1 + term2;
}

double complex ComputeInt(
  gsl_matrix* splinecoeffsAreal,         /*  */
  gsl_matrix* splinecoeffsAimag,         /*  */
  gsl_matrix* quadsplinecoeffsphase)     /*  */
{
  double complex res = 0.;
  double complex resint = 0.;
  /* Number of points - i.e. nb of intervals + 1 */
  /* Assumes that the dimensions match and that the frequency vectors are the same between the different splinecoeffs */
  int nbpts = (int) quadsplinecoeffsphase->size1;

  //clock_t tbegtotal = clock();
  for(int j=0; j<nbpts-1; j++) {
    double eps = gsl_matrix_get(quadsplinecoeffsphase, j+1, 0) - gsl_matrix_get(quadsplinecoeffsphase, j, 0);
    double epspow[4];
    epspow[0] = 1.;
    epspow[1] = eps;
    epspow[2] = eps*eps;
    epspow[3] = eps*epspow[2];
    double p0 = gsl_matrix_get(quadsplinecoeffsphase, j, 1);

    /* Rescale the interval to [0,1] and scale out the constant amplitude term */
    double p1 = gsl_matrix_get(quadsplinecoeffsphase, j, 2) * eps;
    double p2 = gsl_matrix_get(quadsplinecoeffsphase, j, 3) * epspow[2];
    double complex A0 = gsl_matrix_get(splinecoeffsAreal, j, 1) + I*gsl_matrix_get(splinecoeffsAimag, j, 1);
    double complex A0inv = 1./A0;
    double A0abs = cabs(A0); /* Used for accuracy requirement of adaptive Taylor expansions Ia and Ib */
    double complex coeffsA[4] = {0.,0.,0.,0.};
    coeffsA[0] = 1.;
    for(int i=1; i<=3; i++) coeffsA[i] = epspow[i] * A0inv * (gsl_matrix_get(splinecoeffsAreal, j, i+1) + I*gsl_matrix_get(splinecoeffsAimag, j, i+1));

    /* Factor scaled out */
    double complex factor = eps * A0 * cexp(I*p0);

    double absp1 = fabs(p1); double absp2 = fabs(p2);
    if(absp1<p1threshold1 && absp2<p2threshold) {
      resint = factor * ComputeIntCase1a(coeffsA, p1, p2, A0abs);
    }
    else if(p1threshold1<=absp1 && absp1<p1threshold2 && absp2<p2threshold) {
      resint = factor * ComputeIntCase1b(coeffsA, p1, p2, A0abs);
    }
    else if(absp2>=p2p1slopethreshold*absp1 && absp2>=p2threshold) {
      resint = factor * ComputeIntCase2(coeffsA, p1, p2);
    }
    else if(absp2<p2p1slopethreshold*absp1 && absp2>=p2threshold && absp1<p1threshold2) {
      resint = factor * ComputeIntCase3(coeffsA, p1, p2);
    }
    else {
      resint = factor * ComputeIntCase4(coeffsA, p1, p2);
    }
    res += resint;
  }

  return res;
}

max_stars_repo_name: titodalcanton/flare

max_issues_repo_name: titodalcanton/flare

max_forks_repo_name: titodalcanton/flare


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

#ifndef __GSL_MATRIX_COMPLEX_LONG_DOUBLE_H__
#define __GSL_MATRIX_COMPLEX_LONG_DOUBLE_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_complex_long_double.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

typedef struct
{
  size_t size1;
  size_t size2;
  size_t tda;
  long double * data;
  gsl_block_complex_long_double * block;
  int owner;
} gsl_matrix_complex_long_double ;

typedef struct
{
  gsl_matrix_complex_long_double matrix;
} _gsl_matrix_complex_long_double_view;

typedef _gsl_matrix_complex_long_double_view gsl_matrix_complex_long_double_view;

typedef struct
{
  gsl_matrix_complex_long_double matrix;
} _gsl_matrix_complex_long_double_const_view;

typedef const _gsl_matrix_complex_long_double_const_view gsl_matrix_complex_long_double_const_view;

/* Allocation */

GSL_FUN gsl_matrix_complex_long_double *
gsl_matrix_complex_long_double_alloc (const size_t n1, const size_t n2);

GSL_FUN gsl_matrix_complex_long_double *
gsl_matrix_complex_long_double_calloc (const size_t n1, const size_t n2);

GSL_FUN gsl_matrix_complex_long_double *
gsl_matrix_complex_long_double_alloc_from_block (gsl_block_complex_long_double * b,
                                           const size_t offset,
                                           const size_t n1, const size_t n2, const size_t d2);

GSL_FUN gsl_matrix_complex_long_double *
gsl_matrix_complex_long_double_alloc_from_matrix (gsl_matrix_complex_long_double * b,
                                            const size_t k1, const size_t k2,
                                            const size_t n1, const size_t n2);

GSL_FUN gsl_vector_complex_long_double *
gsl_vector_complex_long_double_alloc_row_from_matrix (gsl_matrix_complex_long_double * m,
                                                const size_t i);

GSL_FUN gsl_vector_complex_long_double *
gsl_vector_complex_long_double_alloc_col_from_matrix (gsl_matrix_complex_long_double * m,
                                                const size_t j);

GSL_FUN void gsl_matrix_complex_long_double_free (gsl_matrix_complex_long_double * m);

/* Views */

GSL_FUN _gsl_matrix_complex_long_double_view
gsl_matrix_complex_long_double_submatrix (gsl_matrix_complex_long_double * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

GSL_FUN _gsl_vector_complex_long_double_view
gsl_matrix_complex_long_double_row (gsl_matrix_complex_long_double * m, const size_t i);

GSL_FUN _gsl_vector_complex_long_double_view
gsl_matrix_complex_long_double_column (gsl_matrix_complex_long_double * m, const size_t j);

GSL_FUN _gsl_vector_complex_long_double_view
gsl_matrix_complex_long_double_diagonal (gsl_matrix_complex_long_double * m);

GSL_FUN _gsl_vector_complex_long_double_view
gsl_matrix_complex_long_double_subdiagonal (gsl_matrix_complex_long_double * m, const size_t k);

GSL_FUN _gsl_vector_complex_long_double_view
gsl_matrix_complex_long_double_superdiagonal (gsl_matrix_complex_long_double * m, const size_t k);

GSL_FUN _gsl_vector_complex_long_double_view
gsl_matrix_complex_long_double_subrow (gsl_matrix_complex_long_double * m,
                                 const size_t i, const size_t offset,
                                 const size_t n);

GSL_FUN _gsl_vector_complex_long_double_view
gsl_matrix_complex_long_double_subcolumn (gsl_matrix_complex_long_double * m,
                                    const size_t j, const size_t offset,
                                    const size_t n);

GSL_FUN _gsl_matrix_complex_long_double_view
gsl_matrix_complex_long_double_view_array (long double * base,
                             const size_t n1,
                             const size_t n2);

GSL_FUN _gsl_matrix_complex_long_double_view
gsl_matrix_complex_long_double_view_array_with_tda (long double * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

GSL_FUN _gsl_matrix_complex_long_double_view
gsl_matrix_complex_long_double_view_vector (gsl_vector_complex_long_double * v,
                              const size_t n1,
                              const size_t n2);

GSL_FUN _gsl_matrix_complex_long_double_view
gsl_matrix_complex_long_double_view_vector_with_tda (gsl_vector_complex_long_double * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

GSL_FUN _gsl_matrix_complex_long_double_const_view
gsl_matrix_complex_long_double_const_submatrix (const gsl_matrix_complex_long_double * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

GSL_FUN _gsl_vector_complex_long_double_const_view
gsl_matrix_complex_long_double_const_row (const gsl_matrix_complex_long_double * m,
                            const size_t i);

GSL_FUN _gsl_vector_complex_long_double_const_view
gsl_matrix_complex_long_double_const_column (const gsl_matrix_complex_long_double * m,
                               const size_t j);

GSL_FUN _gsl_vector_complex_long_double_const_view
gsl_matrix_complex_long_double_const_diagonal (const gsl_matrix_complex_long_double * m);

GSL_FUN _gsl_vector_complex_long_double_const_view
gsl_matrix_complex_long_double_const_subdiagonal (const gsl_matrix_complex_long_double * m,
                                    const size_t k);

GSL_FUN _gsl_vector_complex_long_double_const_view
gsl_matrix_complex_long_double_const_superdiagonal (const gsl_matrix_complex_long_double * m,
                                      const size_t k);

GSL_FUN _gsl_vector_complex_long_double_const_view
gsl_matrix_complex_long_double_const_subrow (const gsl_matrix_complex_long_double * m,
                                       const size_t i, const size_t offset,
                                       const size_t n);

GSL_FUN _gsl_vector_complex_long_double_const_view
gsl_matrix_complex_long_double_const_subcolumn (const gsl_matrix_complex_long_double * m,
                                          const size_t j, const size_t offset,
                                          const size_t n);

GSL_FUN _gsl_matrix_complex_long_double_const_view
gsl_matrix_complex_long_double_const_view_array (const long double * base,
                                   const size_t n1,
                                   const size_t n2);

GSL_FUN _gsl_matrix_complex_long_double_const_view
gsl_matrix_complex_long_double_const_view_array_with_tda (const long double * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

GSL_FUN _gsl_matrix_complex_long_double_const_view
gsl_matrix_complex_long_double_const_view_vector (const gsl_vector_complex_long_double * v,
                                    const size_t n1,
                                    const size_t n2);

GSL_FUN _gsl_matrix_complex_long_double_const_view
gsl_matrix_complex_long_double_const_view_vector_with_tda (const gsl_vector_complex_long_double * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

GSL_FUN void gsl_matrix_complex_long_double_set_zero (gsl_matrix_complex_long_double * m);
GSL_FUN void gsl_matrix_complex_long_double_set_identity (gsl_matrix_complex_long_double * m);
GSL_FUN void gsl_matrix_complex_long_double_set_all (gsl_matrix_complex_long_double * m, gsl_complex_long_double x);

GSL_FUN int gsl_matrix_complex_long_double_fread (FILE * stream, gsl_matrix_complex_long_double * m) ;
GSL_FUN int gsl_matrix_complex_long_double_fwrite (FILE * stream, const gsl_matrix_complex_long_double * m) ;
GSL_FUN int gsl_matrix_complex_long_double_fscanf (FILE * stream, gsl_matrix_complex_long_double * m);
GSL_FUN int gsl_matrix_complex_long_double_fprintf (FILE * stream, const gsl_matrix_complex_long_double * m, const char * format);

GSL_FUN int gsl_matrix_complex_long_double_memcpy(gsl_matrix_complex_long_double * dest, const gsl_matrix_complex_long_double * src);
GSL_FUN int gsl_matrix_complex_long_double_swap(gsl_matrix_complex_long_double * m1, gsl_matrix_complex_long_double * m2);

GSL_FUN int gsl_matrix_complex_long_double_swap_rows(gsl_matrix_complex_long_double * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_complex_long_double_swap_columns(gsl_matrix_complex_long_double * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_complex_long_double_swap_rowcol(gsl_matrix_complex_long_double * m, const size_t i, const size_t j);

GSL_FUN int gsl_matrix_complex_long_double_transpose (gsl_matrix_complex_long_double * m);
GSL_FUN int gsl_matrix_complex_long_double_transpose_memcpy (gsl_matrix_complex_long_double * dest, const gsl_matrix_complex_long_double * src);

GSL_FUN int gsl_matrix_complex_long_double_isnull (const gsl_matrix_complex_long_double * m);
GSL_FUN int gsl_matrix_complex_long_double_ispos (const gsl_matrix_complex_long_double * m);
GSL_FUN int gsl_matrix_complex_long_double_isneg (const gsl_matrix_complex_long_double * m);
GSL_FUN int gsl_matrix_complex_long_double_isnonneg (const gsl_matrix_complex_long_double * m);

GSL_FUN int gsl_matrix_complex_long_double_add (gsl_matrix_complex_long_double * a, const gsl_matrix_complex_long_double * b);
GSL_FUN int gsl_matrix_complex_long_double_sub (gsl_matrix_complex_long_double * a, const gsl_matrix_complex_long_double * b);
GSL_FUN int gsl_matrix_complex_long_double_mul_elements (gsl_matrix_complex_long_double * a, const gsl_matrix_complex_long_double * b);
GSL_FUN int gsl_matrix_complex_long_double_div_elements (gsl_matrix_complex_long_double * a, const gsl_matrix_complex_long_double * b);
GSL_FUN int gsl_matrix_complex_long_double_scale (gsl_matrix_complex_long_double * a, const gsl_complex_long_double x);
GSL_FUN int gsl_matrix_complex_long_double_add_constant (gsl_matrix_complex_long_double * a, const gsl_complex_long_double x);
GSL_FUN int gsl_matrix_complex_long_double_add_diagonal (gsl_matrix_complex_long_double * a, const gsl_complex_long_double x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
GSL_FUN int gsl_matrix_complex_long_double_get_row(gsl_vector_complex_long_double * v, const gsl_matrix_complex_long_double * m, const size_t i);
GSL_FUN int gsl_matrix_complex_long_double_get_col(gsl_vector_complex_long_double * v, const gsl_matrix_complex_long_double * m, const size_t j);
GSL_FUN int gsl_matrix_complex_long_double_set_row(gsl_matrix_complex_long_double * m, const size_t i, const gsl_vector_complex_long_double * v);
GSL_FUN int gsl_matrix_complex_long_double_set_col(gsl_matrix_complex_long_double * m, const size_t j, const gsl_vector_complex_long_double * v);
/***********************************************************************/

/* inline functions if you are using GCC */

GSL_FUN INLINE_DECL gsl_complex_long_double gsl_matrix_complex_long_double_get(const gsl_matrix_complex_long_double * m, const size_t i, const size_t j);
GSL_FUN INLINE_DECL void gsl_matrix_complex_long_double_set(gsl_matrix_complex_long_double * m, const size_t i, const size_t j, const gsl_complex_long_double x);

GSL_FUN INLINE_DECL gsl_complex_long_double * gsl_matrix_complex_long_double_ptr(gsl_matrix_complex_long_double * m, const size_t i, const size_t j);
GSL_FUN INLINE_DECL const gsl_complex_long_double * gsl_matrix_complex_long_double_const_ptr(const gsl_matrix_complex_long_double * m, const size_t i, const size_t j);

#ifdef HAVE_INLINE

INLINE_FUN
gsl_complex_long_double
gsl_matrix_complex_long_double_get(const gsl_matrix_complex_long_double * m,
                     const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      gsl_complex_long_double zero = {{0,0}};

      if (i >= m->size1)
        {
          GSL_ERROR_VAL("first index out of range", GSL_EINVAL, zero) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VAL("second index out of range", GSL_EINVAL, zero) ;
        }
    }
#endif
  return *(gsl_complex_long_double *)(m->data + 2*(i * m->tda + j)) ;
}

INLINE_FUN
void
gsl_matrix_complex_long_double_set(gsl_matrix_complex_long_double * m,
                     const size_t i, const size_t j, const gsl_complex_long_double x)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_VOID("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VOID("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  *(gsl_complex_long_double *)(m->data + 2*(i * m->tda + j)) = x ;
}

INLINE_FUN
gsl_complex_long_double *
gsl_matrix_complex_long_double_ptr(gsl_matrix_complex_long_double * m,
                             const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (gsl_complex_long_double *)(m->data + 2*(i * m->tda + j)) ;
}

INLINE_FUN
const gsl_complex_long_double *
gsl_matrix_complex_long_double_const_ptr(const gsl_matrix_complex_long_double * m,
                                   const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (const gsl_complex_long_double *)(m->data + 2*(i * m->tda + j)) ;
}

#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_MATRIX_COMPLEX_LONG_DOUBLE_H__ */

ext: h

max_stars_repo_name: berkus/music-cs

max_issues_repo_name: berkus/music-cs

max_forks_repo_name: berkus/music-cs