/********************************************************************************
 * Project: Exact Bayesian Regression of Piecewise Constant Functions           *
 * Module : pcreg.cpp                                                           *
 * Author : MH                                                                  *
 * Content:                                                                     *
 * (c) 2001 by Marcus Hutter                                                    *
 ********************************************************************************/

#include <float.h>   // DBL_MAX
#include <limits.h>   // INT_MAX
#include <windows.h> // BOOL
#include <stdlib.h>
#include <stdio.h> // size_t
#include <string.h> // size_t
#include <math.h>

/********************************************************************************/
/*                        D e s c r i p t i o n                                 */
/********************************************************************************

This module implements exact and efficient Bayesian regression for
piecewise constant functions of unknown segment number, boundary
location, and levels. It works for any noise and segment level
prior, e.g. Cauchy which can handle outliers. Simple
but good estimates for the in-segment variance are computed.
A Bayesian regression curve as a better way of smoothing data
without blurring boundaries is also computed. The Bayesian approach also allows
determination of the evidence, break probabilities
and error estimates, useful for model selection and significance
and robustness studies. The module contains further synthetic data generation
and loading of data from files. Output is to a file.

This model can be downloaded from website
http://www.idsia.ch/~marcus/ai/pcreg.htm
which also contains a detaild description of its working.

Functions:
----------
- LogAdd()              Overflow-safe computation of log(e^x+e^y)
- lnGamma()             log(Gamma(n)) - by pretabulation
- lnBinom()             log of Binomial Coefficient (n k)
- GaussDistr()          returns Gauss(x|mean=n,variance=v)
- CauchyDistr()         returns Cauchy(x|"mean"=m,"variance"=v)
- SampleGauss()         Sample from Gaussian Distribution with mean=0 and variance=1
- SampleCauchy()        Sample from Cauchy Distribution p(x)=1/(1+x^2)
- ReadDataCol()         Read data column from file 
- ReadKeyDataCol()      Read data column from file by key
- CreatePCData()        Create synthetic data or read data from file
- EstGlobParam()        Estimate global parameters mglob, varglob, varseg from data y
- CrSegDist1Gauss()     Create Gaussian in-segment evidence, mean, and variance
- CrSegDist1Cauchy()    Create Cauchy in-segment evidence, mean, and variance
- CrSegDist1()          Create in-segment evidence, mean, and variance based on rmode
- PCEvidence()          Determine Left and Right Evidences lL and lR from lA0 
- RegMoms()             Bayesian regression curve (moments) to S
- PCReg()               Determine PC-regression and regression curve
- PCEvdKml()            Computes evidence and kml for data y (fixed varglob, ...)
- PCLikelihood()        Computes data log-likelihood ll and expected likelihood ElogP with variance VlogP
- prPCRResults()        Print regression results to dev
- PCRegPrint()          Perform regression and print results to dev and dev2 
- PCRegGenomes()        Regression for 23 Chromosomes of 31 Patients
/- EstPrintDensity()    Estimate and print distribution of data
- EvdFctSeg()           Plot Evidence and kml as function of segment variance varseg
- EvdFctGlob()          Plot Evidence and kml as function of segment variance varglob
- Exp4Paper()           Batch Regression for all Examples in Paper

Notation:
---------
- double y[i]      = data item i (i=1..n)
- int    n         = data size
- double lA0[i][j] = log-evidence of single segment i+1..j (0<=i<j<=n)
- double A10[i][j] = mean of single segment i+1..j
- double A20[i][j] = 2nd moment of single segment i+1..j
- double lL[k][j]  = log-evidence of k segments from 1..j
- double lR[k][i]  =  log-evidence of k segments from i+1..n
- int    kmax      = upper bound on number of segments (1<=kmax<=n)
- double lC[k]     = log P(k|y) = log-probability of k segments
- double lE        = log P(y) = log-evidence of data y
- int    kml       = maximum likelihood (ML) estimate of k
- double B[i]      = Probability of break at i (i=0..n)
- int    J[k]      = end of segment k (0<=k<=kml)
- double M1[k]     = (mean) level of segment k
- double Mt[i]     = piecwise constant regression curve
- double S[1..n]   = Bayesian regression curve
- double V[i]      = estimated variance of curve at point i

Adapt global path "../../../Temp/Regression/" for I/O

*/

/********************************************************************************/
/*                   E l e m e n t a r y   F u n c t i o n                      */
/********************************************************************************/

#define rrand() (((double)(rand()+(rand()<<15)+1))/((1<<30)+1)) /* random number in ]0..1[          */
#define random2(n)      ((int)(((long)rand()*(n))>>15)) /* integer random number in 0..n-1          */
#define PI   3.141592653589793238462
FILE *dev=0;               // standard output device
FILE *dev2=0;              // tabular material in extra file (or same dev2=dev)
const int maxnprint=10000; // length at which printf tables will be cut
#ifdef _DEBUG
const int check=1;         // check but slow
#else
const int check=0;         // no check but fast
#endif

/*------------------------------------------------------------------------------*
  Description:  Error handler
 *-MH---------------------------------------------------------------------------*/
int panicf(char *s) 
/*------------------------------------------------------------------------------*/
{ 
  printf(s);  // set a breakpoint here
  return 0;
}  // End of panicf()

/*------------------------------------------------------------------------------*
  Description:  Allocate memory for matrix A[0..N-1][0..M-1] with components 
                of type typ and size typl. A is a vector of pointers to vectors.
  Parameters:   N    - number of columns.
                M    - number of rows.
                typ  - type of vector elements (e.g. int)
                typl - sizeof(typ).
  Return Value: pointer to allocated matrix or NULL if allocation error.
  Remark1:      Use macro Matrix().
  Remark2:      A is a vector of pointers to vectors. This allows to use A as
                parameter because number of columns need not to be known.
  Remark3:      free space with standard free() function.
 *-MH---------------------------------------------------------------------------*/
void*   XMAlloc(int N, int M, size_t typl)
/*------------------------------------------------------------------------------*/
{
  int   i;
  size_t VS=M*typl;                     /* size of row-vector                   */
  char **lA=(char**)malloc(N*sizeof(void*) + N*VS); // void *lA=...
  if (lA==NULL)
  {  panicf("allocation failure in 'Matrix'");
     return(NULL);
  }
  lA[0]=(char*)(lA+N);               //((void**)lA)[0]=((void**)lA)+N;
  for(i=1;i<N;i++) lA[i]=lA[i-1]+VS; //((void**)lA)[i]=((BYTE**)lA)[i-1]+VS;
  return((void*)lA);
}

#define CrMatrix(N,M,typ) ((typ**)XMAlloc(N,M,sizeof(typ)))

#define CrVector(N,typ)   ((typ*)malloc((N)*sizeof(typ)))

/*------------------------------------------------------------------------------*
  Description:  Overflow-safe computation of log(e^x+e^y)
 *-MH---------------------------------------------------------------------------*/
double LogAdd(double x, double y)
/*------------------------------------------------------------------------------*/
{
  if (x>y+40) return x;                /* just speedup, exact for double        */
  if (y>x+40) return y;                /* just speedup, exact for double        */
  if (x<y) return y+log(1+exp(x-y));
  else     return x+log(1+exp(y-x));
} // End of LogAdd()

/*------------------------------------------------------------------------------*
  Description:  log(Gamma(n)) - by pretabulation
  Remark:       Currently only implemented for integer n
 *-MH---------------------------------------------------------------------------*/
double lnGamma(int n)
/*------------------------------------------------------------------------------*/
{
  static double lnGtab[2002] = { -1 };
  if (n>2001) panicf("lnGamma() not implemented for n>2001\n");

  if (lnGtab[0]<0)
  { lnGtab[0]=DBL_MAX; lnGtab[1]=0;
    for(int i=1;i<2001;i++) lnGtab[i+1]=lnGtab[i]+log(i);
  }
  return lnGtab[n];
} // End of lnGamma()

/*------------------------------------------------------------------------------*
  Description:  log of Binomial Coefficient (n k)
  Return Value: log(n!/k!/(n-k)!)
 *-MH---------------------------------------------------------------------------*/
double lnBinom(int n, int k)
/*------------------------------------------------------------------------------*/
{
  return lnGamma(n+1)-lnGamma(k+1)-lnGamma(n-k+1);
} // End of lnBinom()

/*------------------------------------------------------------------------------*
  Description:  returns Gauss(x|mean=n,variance=v)
 *-MH---------------------------------------------------------------------------*/
double GaussDistr(double x, double m, double v)
/*------------------------------------------------------------------------------*/
{
  return exp(-(x-m)*(x-m)/(2*v))/sqrt(2*PI*v);
} // End of GaussDistr()

/*------------------------------------------------------------------------------*
  Description:  returns Cauchy(x|"mean"=m,"variance"=v)
 *-MH---------------------------------------------------------------------------*/
double CauchyDistr(double x, double m, double v)
/*------------------------------------------------------------------------------*/
{
  return sqrt(v)/(v+(x-m)*(x-m))/PI;
} // End of CauchyDistr()

/*------------------------------------------------------------------------------*
  Description:  Sample from Gaussian Distribution with mean=0 and variance=1
 *-MH---------------------------------------------------------------------------*/
double SampleGauss()
/*------------------------------------------------------------------------------*/
{
  static int iset=0;
  static double gset;
  double fac,rsq,v1,v2;

  if  (iset==0) 
  { do 
    { v1=2*rrand()-1;
      v2=2*rrand()-1;
      rsq=v1*v1+v2*v2;
    } while (rsq>=1||rsq==0);
    fac=sqrt(-2*log(rsq)/rsq);
    gset=v1*fac;
    iset=1;
    return v2*fac;
  }else                                                             
  { iset=0;
    return gset;
  }
} // SampleGauss()

/*------------------------------------------------------------------------------*
  Description:  Sample from Cauchy Distribution p(x)=1/(1+x^2)
 *-MH---------------------------------------------------------------------------*/
double SampleCauchy()
/*------------------------------------------------------------------------------*/
{
  return tan(PI*(rrand()-0.5));
} // SampleCauchy()


/********************************************************************************/
/*                          D a t a   G e n e r a t i o n                       */
/********************************************************************************/

double varseg;                         // standard noise variance in segment (sigma^2)
double varglob;                        // global segment level variance (rho^2)
double mglob;                          // global segment level mean (nu)
int    nmax  = 1000;                   // maximal data size / or exact data size for synthetic
int    kmax  = min(nmax,100);          // maximal number of segments
int    rmode = 1;                      // regression mode (1=Gauss=fast, 2=Cauchy=slow)
int    dmode = 4;                      // choice of synthetic data dName[data]
const char *rName[] = { "---", "Gauss", "Cauchy", "Non-Parametric", "other" };
const char *dName[] = { "Test", "Uniform", "Gauss3", "Cauchy3", "GeneExpression", "Gauss3[_~--]", "Cauchy3[_~--]" };
double vsfudge = 1.0;                  // fudge multiplier of estimate of varseg 
double mixk = 0;                       // 0 = B,S,V for kml, 1 = B,S,V mixed over k
double unbiased = 1;                   // 0 = exact S,V, 1 = unbiased S,V (Gauss only)

int    datacol=5;                      // data column to analyse (default=3)
int    datakey=9;                      // only data with key (default=1)
int    gcnt=1;                         // global data item counter (only used for print)

/*------------------------------------------------------------------------------*
  Description:  Read data column from file 
  Parameters:   Read column col from row row1st to row1st+n-1 
                from file dName to array y[1..n] 
  Return Value: maximal number of n successful read
  Remark:       
 *-MH---------------------------------------------------------------------------*/
int ReadDataCol(char *dName, double y[], int n, int row1st, int col)
/*------------------------------------------------------------------------------*/
{
  char   c;
  FILE *devr=fopen(dName,"r");
  if (!devr) panicf("Can't open data file for read\n");
  for(int t=1;t<row1st;t++) do c=fgetc(devr); while(c!='\n'&& c!=EOF);
  for(t=1;t<=n;t++)
  {
    for(int i=1;i<col;i++) fscanf(devr,"%*s");
    if (fscanf(devr,"%lf",&y[t])==EOF) { n=t-1; break; }
    do c=fgetc(devr); while(c!='\n'&& c!=EOF);
  }
  fclose(devr);
  return n;
} // End of ReadDataCol()

/*------------------------------------------------------------------------------*
  Description:  Read data column from file by key
  Parameters:   Read column col from rows which contain key in column keycol
                from file dName to array y[1..n] 
  Return Value: number of n successful read 
  Remark:       Works only for keycol<col
 *-MH---------------------------------------------------------------------------*/
int ReadKeyDataCol(char *dName, double y[], int n, int col, int key, int keycol)
/*------------------------------------------------------------------------------*/
{
  char   c,temp[99];
  int    i,t,rkey;
  FILE *devr=fopen(dName,"r");
  if (!devr) panicf("Can't open data file for read\n");
  for(t=1;;t++)
  {
    for(i=1;i<keycol;i++) fscanf(devr,"%s",temp);
    fscanf(devr,"%d",&rkey);
    if (rkey!=key) { do c=fgetc(devr); while(c!='\n'&& c!=EOF); }
    else break;
  }
  gcnt=t; // optional
  for(t=1;t<=n;t++)
  {  
    for(i=keycol+1;i<col;i++) fscanf(devr,"%s",temp);
    if (fscanf(devr,"%lf",&y[t])==EOF) break;
    do c=fgetc(devr); while(c!='\n'&& c!=EOF);

    for(i=1;i<keycol;i++) fscanf(devr,"%*s");
    fscanf(devr,"%d",&rkey);
    if (rkey!=key) { t++; break; }
  }
  fclose(devr);
  return n=t-1;
} // End of ReadKeyDataCol()

/*------------------------------------------------------------------------------*
  Description:  Create synthetic data or read data from file
  Parameters:   in: n = (upper bound on) data size 
                out: array y[0..n]
  Return Value: data size (possibly <n)
  Remark:       For synthetic functions varseg is known and may be used for validation
 *-MH---------------------------------------------------------------------------*/
int CreatePCData(double y[], int n, int dmode)
/*------------------------------------------------------------------------------*/
{
  double stdsegm;                      /* true std-deviation of segment in model*/
  int i;
  switch(dmode)
  {
    case 0: // Test
      stdsegm=0.3;
      for(i=1;i<=n;i++) y[i]=stdsegm*SampleGauss();
      break;
    case 1: // uniform [0,1]
      for(i=1;i<=n;i++) y[i]=rrand();
//      varseg=1/12.0;
      break;
    case 2: // Three equi-Gaussian segments (within model-assumption)
      stdsegm=1;
      for(i=1;i<=n;i++) y[i]=stdsegm*SampleGauss();
      for(i=1;i<=n/3;i++)   y[i]+=1;
      for(i=2*n/3+1;i<=n;i++) y[i]-=1;
//      varseg=stdsegm*stdsegm;
      break;
    case 3: // Three equi-Cauchy segments (to model outliers)
      stdsegm=0.1;
      for(i=1;i<=n;i++) y[i]=stdsegm*SampleCauchy();
      for(i=1;i<=n/3;i++)   y[i]+=1;
      for(i=2*n/3+1;i<=n;i++) y[i]-=1;
//      varseg=stdsegm*stdsegm; // 0.25*0.25; if Gauss is used for fit
      break;
    case 4: // IOSI mm-data (datacol,keycol=2,datakey=gene#)
//      ReadDataColAv("mm-copynumber-normalized.dat",y,n,2,2,1);
      n=ReadKeyDataCol("../../../Temp/Regression/mm-copynumber-normalized.dat",y,n,datacol,datakey,2);
      for(i=1;i<=n;i++) y[i]=log(y[i]);
//      varseg=0.1; // Column lC (now D X05.090_SPA_CN) is one perfect Gauss(mean=0.69,var=0.087) segment
      break;
    case 5: // Three Gaussian segments [_~--]
      for(i=1;i<=n;i++) y[i]=mglob+sqrt(varseg)*SampleGauss();
      for(i=1;i<=n/4;i++) y[i]-=sqrt(varglob);
      for(i=n/4+1;i<=n/2;i++) y[i]+=sqrt(varglob);;
      break;
    case 6: // Three Cauchy segments [_~--]
      for(i=1;i<=n;i++) y[i]=mglob+sqrt(varseg)*SampleCauchy();
      for(i=1;i<=n/4;i++) y[i]-=sqrt(varglob);
      for(i=n/4+1;i<=n/2;i++) y[i]+=sqrt(varglob);;
      break;
    default:
      panicf("Wrong function (dmode)\n");
  }
  return n;
} // End of CreatePCData()

/********************************************************************************/
/*                         P C - R e g r e s s i o n                            */
/********************************************************************************/

/*------------------------------------------------------------------------------*
  Description:  Compare two double's *a and *b (needed for qsort)
 *-MH---------------------------------------------------------------------------*/
int DoubleComp(const void *a, const void *b)
/*------------------------------------------------------------------------------*/
{
  double d=(*(double*)a) - (*(double*)b);
  return d==0?0:(d<0?-1:1);
} // End of DoubleComp()

/*------------------------------------------------------------------------------*
  Description:  Estimate global parameters mglob, varglob, varseg from data y
  Remark:       For synthetic functions varseg is known and may be used instead for validation
 *-MH---------------------------------------------------------------------------*/
void EstGlobParam(double y[], int n)
/*------------------------------------------------------------------------------*/
{
  if (rmode==1)
  { double m=0,s=0,l=0;
    for(int i=1;i<=n;i++) { m+=y[i]; s+=y[i]*y[i]; }
    for(i=1;i<n;i++) l+=(y[i]-y[i+1])*(y[i]-y[i+1]);
    mglob=m/n;
    varglob=s/n-(m/n)*(m/n);
    varseg=vsfudge*l/2/(n-1);
    if (varglob-varseg<varglob/10) printf("varglob<<varseg ==> strange data in EstGlobParam()\n");
//    varglob=max(varglob/10,varglob-varseg); //  more accurate (but bad and could be negative)
  }else
  { const double alpha=1;              // 1=Cauchy, 0.67448975=Gauss
    const double beta=2;               // 2=Cauchy, 0.67448975*sqrt(2)=Gauss
    double *ys=CrVector(n+1,double);
    memcpy(ys,y,(n+1)*sizeof(double));
    qsort(ys+1,n,sizeof(double),DoubleComp);
    mglob=ys[n/2+1];
    varglob=(ys[(3*n)/4+1]-ys[n/4+1])/2/alpha;
    for(int t=1;t<n;t++) ys[t]=y[t+1]-y[t];
    qsort(ys+1,n-1,sizeof(double),DoubleComp);
    varseg=(ys[(3*n)/4+1]-ys[n/4+1])/2/beta;
    if (varglob-varseg<varglob/10) printf("varglob<<varseg ==> strange data in EstGlobParam()\n");
//    varglob=max(varglob/10,varglob-varseg); //  more accurate (but bad and could be negative)
    // square after shift, since Cauchy is additive in std, not var
    varglob=varglob*varglob;
    varseg=varseg*varseg;
  }
} // End of EstGlobParam()

/*------------------------------------------------------------------------------*
  Description:  Create Gaussian in-segment evidence, mean, and variance
  Parameters:   in: y, n, mglob, varglob, varseg
                out: lA0, A10 (may be NULL), A20 (may be NULL if also A10 is)
  Remark:       Verified against CrSegDist1old() and CrSegDistCauchy()
 *-MH---------------------------------------------------------------------------*/
void CrSegDist1Gauss(double y[], int n, double **lA0, double **A10, double **A20)
/*------------------------------------------------------------------------------*/
{
  if (n>999) printf("CrSegDist1Gauss(): Determine 1st Order Distance-Matrix & Model Dimension [%d*]:\n",n);
  for(int i=0;i<=n;i++)
  { double ysum=0,y2sum=0;
    if (n>999) printf("*");                      
    for(int j=0;j<=n;j++)
    {
      int d=j-i;
      if (d<=0) { lA0[i][j]=-DBL_MAX; if (A10) A10[i][j]=0; if (A20) A20[i][j]=0; continue; }
      ysum+=y[j]-mglob; y2sum+=(y[j]-mglob)*(y[j]-mglob);
      lA0[i][j]=-0.5*d*log(2*PI*varseg)-0.5*log(1+d*varglob/varseg)
//                +0.5/varseg*(ysum*ysum/(d+varseg/varglob)-y2sum); // equivalent
                -0.5/varseg*(varglob*(d*y2sum-ysum*ysum)+varseg*y2sum)/(d*varglob+varseg); // equivalent
      if (fabs(ysum*ysum/(d+varseg/varglob)-y2sum + (varglob*(d*y2sum-ysum*ysum)+varseg*y2sum)/(d*varglob+varseg))>1e-9) printf("Incons in Cr\n");
      if (A10&&!unbiased) A10[i][j]=(varglob*(ysum+d*mglob)+varseg*mglob)/(d*varglob+varseg); // exact
      if (A20&&!unbiased) A20[i][j]=A10[i][j]*A10[i][j]+1/(d/varseg+1/varglob);               // exact
      if (A10&&unbiased) A10[i][j]=mglob+ysum/d;                 // unbiased
      if (A20&&unbiased) A20[i][j]=A10[i][j]*A10[i][j]+varseg/d; // unbiased
    }
  } 
  if (n>999) printf("\n");
} // End of CrSegDist1Gauss()

/*------------------------------------------------------------------------------*
  Description:  Create Cauchy in-segment evidence, mean, and variance
  Parameters:   in: y, n, mglob, varglob, varseg
                out: lA0, A10 (may be NULL), A20 (may be NULL if also A10 is)
  Remark:       Verified against CrSegDist1Gauss()
                slow: time O(n^2*integral)
                Simply replace CauchyDistr by other distributions when needed
 *-MH---------------------------------------------------------------------------*/
void CrSegDist1Cauchy(double y[], int n, double **lA0, double **A10, double **A20)
/*------------------------------------------------------------------------------*/
{
  double irange=50*sqrt(varglob);
  double istep=sqrt(varseg)/10;
  int    inum=(int)(irange/istep);
  double istart=mglob-irange/2;
  double *lPymu=CrVector(inum,double);
//  double lPymu[99];

  for(int i=0;i<=n;i++)
  { 
    printf("*");
    for(int l=0;l<inum;l++) lPymu[l]=log(CauchyDistr(istart+l*istep,mglob,varglob)*istep);
    for(int j=0;j<=n;j++)
    {
      int d=j-i;
      double s=-DBL_MAX; 
      if (d<=0) 
      { lA0[i][j]=-DBL_MAX; 
        if (A10) A10[i][j]=0; 
        if (A20) A20[i][j]=0; 
        if (check&&d==0) for(l=0;l<inum;l++) s=LogAdd(s,lPymu[l]);
        if (check&&d==0&&fabs(s)>0.1) printf("Integration error larger than %lf in CrSegDist1Cauchy()\n",fabs(s));
        continue; 
      }
      for(l=0;l<inum;l++) //slow?
      { lPymu[l]+=log(CauchyDistr(y[j],istart+l*istep,varseg));
        s=LogAdd(s,lPymu[l]);
      }
      lA0[i][j]=s;
      if (A10) { A10[i][j]=0; for(l=0;l<inum;l++) A10[i][j]+=(istart+l*istep)*exp(lPymu[l]-s); }
      if (A20) { A20[i][j]=0; for(l=0;l<inum;l++) A20[i][j]+=(istart+l*istep)*(istart+l*istep)*exp(lPymu[l]-s); }
    }
  } 
  printf("\n");  
  free(lPymu);
} // End of CrSegDist1Cauchy()

/*------------------------------------------------------------------------------*
  Description:  Create in-segment evidence, mean, and variance based on rmode
  Parameters:   in: y, n, mglob, varglob, varseg
                out: lA0, A10 (may be NULL), A20 (may be NULL if also A10 is)
 *-MH---------------------------------------------------------------------------*/
void CrSegDist1(double y[], int n, double **lA0, double **A10, double **A20)
/*------------------------------------------------------------------------------*/
{
  if (rmode==1) { CrSegDist1Gauss(y,n,lA0,A10,A20); return; }
  if (rmode==2) { CrSegDist1Cauchy(y,n,lA0,A10,A20); return; }
} // End of CrSegDist1()

/*------------------------------------------------------------------------------*
  Description:  Determine Left and Right Evidences lL and lR from lA0 
                up to at most kmax segments
  Return Value: true kmax
  Remark:       Verified for n=3, if called multiply for same data
                slow: time O(kmax*n^2)
                Has some heuristic to stop before kmax if evidences are too small
 *-MH---------------------------------------------------------------------------*/
int PCEvidence(double **lA0, int n, double **lL, double **lR, int kmax)
/*------------------------------------------------------------------------------*/
{
  if (n>100) printf("PCEvidence(): Determine Order 1..kmax Distance Matrix [%d*]:\n",kmax);
  for(int j=0;j<=n;j++) lL[0][j]=lR[0][j]=-DBL_MAX;
  lL[0][0]=lR[0][n]=0;
  int kmx=kmax; double c,cmax=-DBL_MAX;
  for(int k=0;k<kmx;k++)
  {
    if (n>100) printf("*");
    for(j=0;j<=n;j++)
    {
      lL[k+1][j]=lR[k+1][j]=-DBL_MAX;
      for(int h=k;h<=j-1;h++) lL[k+1][j]=LogAdd(lL[k+1][j],lL[k][h]+lA0[h][j]);
      for(h=j+1;h<=n-k;h++)   lR[k+1][j]=LogAdd(lR[k+1][j],lA0[j][h]+lR[k][h]);
    }
    if (fabs(lL[k+1][n]-lR[k+1][0])>1e-9) panicf("Inconsistency in PCEvidence()\n");
    c=lL[k+1][n]-lnBinom(n-1,k);          // 
    if (c>cmax) cmax=c;                   // lowers kmax if possible for speedup
    if (c>cmax-5) kmx=min(kmax,3*k/2+10); //
  }
  if (n>100) printf("\n");
  return kmx;  
} // End of PCEvidence()

/*------------------------------------------------------------------------------*
  Description:  Bayesian regression curve (moments) to S
  Parameters:   Ar0=NULL: normalization/consistency check 
                Ar0=A10: mean is computed
                Ar0=A20: 2nd moment is computed 
  Remark:       time O(kml*n^2)
 *-MH---------------------------------------------------------------------------*/
void RegMoms(int n, int kml, int kmax, double lE, double **lL, double **lR, double **lD, double **lA0, double **Ar0, double S[])
/*------------------------------------------------------------------------------*/
{
  int i,j,h,k;
  double **F=CrMatrix(n+1,n+1,double);
  for(i=0;i<=n;i++) 
  { if (n>100) printf("*"); 
    for(j=0;j<=n;j++)
    {
      F[i][j]=0;
      if (!mixk) for(k=1;k<=kml;k++) F[i][j]+=(Ar0?Ar0[i][j]:1)*exp(lL[k-1][i]+lA0[i][j]+lR[kml-k][j]-lL[kml][n]);
      else      for(k=1;k<=kmax;k++) F[i][j]+=(Ar0?Ar0[i][j]:1)*exp(lL[k-1][i]+lA0[i][j]+lD[k][j]-lE)/kmax;
      if (i>=j&&fabs(F[i][j])>1e-9) panicf("F-Inconsistency in RegMoms()\n");
  } }
  S[0]=0;
  for(h=0;h<n;h++) 
  {
    S[h+1]=S[h]; 
    for(j=h;j<=n;j++) S[h+1]+=F[h][j];
    for(i=0;i<h;i++)  S[h+1]-=F[i][h];
    // verify correctness
    if (!Ar0&&fabs(S[h+1]-1)>1e-9) printf("S0-Inconsistency in RegMoms() 1!=S[%d]=%lf\n",h,S[h]);
    if (check&&n<110) // check against simple and slow O(n^3)
    { printf("*"); double Sh=0; for(i=0;i<h;i++) for(j=h;j<=n;j++) Sh+=F[i][j];
      if (fabs(S[h]-Sh)>1e-9) printf("S-Inconsistency in RegMoms() S[%d]=%lf\n",h,S[h]);
    } 
  }
  if (n>100) printf("\n"); 
  free(F);
} // End of RegMoms()

/*------------------------------------------------------------------------------*
  Description:  Determine PC-regression and regression curve
  Input:        lA0, A10|0, A20|0, lL, lR, n
  Output:       lC, lE, kml, B|0, J|0, S if A10|0, V if A20|0
  Return Value: kml
 *-MH---------------------------------------------------------------------------*/
int PCReg(double **lA0, double **A10, double **A20, double **lL, double **lR, int n, int kmax, 
                 double lC[], double &lE, int &kml, double B[], int J[], double S[], double V[])
/*------------------------------------------------------------------------------*/
{
  int h,k;
  lE=-DBL_MAX;
  kml=1;                // kml based on ML number of segments k

  // Determine evidence lE and lC[k]=p(k|y) and MAP kml=argmax lC[k]
  if (n>100) printf("PCReg(): Determine Regression [%d*]:\n",2*n);
  for(k=1;k<=kmax;k++)
  {
    lC[k]=lL[k][n]-lnBinom(n-1,k-1);
    lE=LogAdd(lE,lC[k]);
    if (lC[k]>lC[kml]) kml=k;
  }
  for(k=1;k<=kmax;k++) lC[k]-=lE;
  lE-=log(kmax);
  // lD matrix
  double **lD=CrMatrix(kmax+1,n+1,double);
  if (mixk)
  { 
    for(int p=0;p<=kmax;p++) for(h=0;h<=n;h++)
    { lD[p][h]=-DBL_MAX;
      for(k=p;k<=kmax;k++) lD[p][h]=LogAdd(lD[p][h],lR[k-p][h]-lnBinom(n-1,k-1));
    }
    // Boundary probabilities (mixture over k)
    if (B)
    { double norm=0;
      for(h=0;h<=n;h++)
      { B[h]=0;
        for(k=1;k<=kmax;k++) 
           B[h]+=exp(lL[k][h]+lD[k][h]-lE)/kmax;
        norm+=B[h];
      }  
      for(k=0;k<=kmax;k++) norm-=k*exp(lC[k]);
      if (fabs(norm)>1e-9) panicf("BreakPoint Inconsistency in PCReg()\n");
      B[0]=1;
    }
  }
  // Regression boundaries (for kml)
  double Bkh;
  if (B&&!mixk) for(h=0;h<=n;h++) B[h]=0; // Boundary probabilities for kml
  for(k=0;k<=kml;k++) // for all segment-boundaries 1..kml-1
  {
    double norm=0;
    double Bkhmax=0;
    for(h=0;h<=n;h++)
    {
      Bkh=exp(lL[k][h]+lR[kml-k][h]-lL[kml][n]);
      if (Bkh>Bkhmax) { Bkhmax=Bkh; if (J) J[k]=h; }
      norm+=Bkh;
      if (B&&!mixk) B[h]+=Bkh; // Boundary probabilities for kml
    }
    if (fabs(norm-1)>1e-9) panicf("Normalization Error in PCReg()\n");
  }
  if (J&&(J[0]!=0||J[kml]!=n)) panicf("Inconsistency in PCReg()\n");
  // Bayesian regression curve and variance curve
  if (S&&check) RegMoms(n,kml,kmax,lE,lL,lR,lD,lA0,0,S); // strong consistency check
  if (S) RegMoms(n,kml,kmax,lE,lL,lR,lD,lA0,A10,S);
  if (V) RegMoms(n,kml,kmax,lE,lL,lR,lD,lA0,A20,V);
  if (V) for(h=1;h<=n;h++) V[h]-=S[h]*S[h];              
  free(lD);
  return kml;  
} // End of PCReg()

/*------------------------------------------------------------------------------*
  Description:  Computes evidence and kml for data y (fixed varglob, ...)
 *-MH---------------------------------------------------------------------------*/
int PCEvdKml(double y[], int n, double &lE, int &kml)
/*------------------------------------------------------------------------------*/
{
  int kmax=min(n,100);
  double **lA0=CrMatrix(n+1,n+1,double);
  CrSegDist1(y,n,lA0,0,0);
  double **lL=CrMatrix(kmax+1,n+1,double);
  double **lR=CrMatrix(kmax+1,n+1,double);
  kmax=PCEvidence(lA0,n,lL,lR,kmax); // can make kmax smaller if appropriate
  double *lC=CrVector(kmax+1,double);
  PCReg(lA0,0,0,lL,lR,n,kmax,lC,lE,kml,0,0,0,0);
  free(lC); free(lR); free(lL); free(lA0); 
  return 1;
} // End of PCEvdKml()

/*------------------------------------------------------------------------------*
  Description: Computes data log-likelihood ll and expected likelihood ElogP with variance VlogP
 *-MH---------------------------------------------------------------------------*/
int PCLikelihood(double y[], int n, double Mt[], double varseg, double &ll, double &ElogP, double &VlogP)
/*------------------------------------------------------------------------------*/
{
  ll=0;
  for(int t=1;t<=n;t++)
  {
    if (rmode==1) ll+=log(GaussDistr(y[t],Mt[t],varseg));
    else          ll+=log(CauchyDistr(y[t],Mt[t],varseg));
  }
  if (rmode==1) { ElogP=-n*(1+log(2*PI*varseg))/2; VlogP=0.5*n; }
  else          { ElogP=-n*log(4*PI*sqrt(varseg)); VlogP=n*PI*PI/3; }
  return 1;
} // End of PCLikelihood()

/********************************************************************************/
/*                P r i n t   R e g r e s s i o n   R e s u l t s               */
/********************************************************************************/

/*------------------------------------------------------------------------------*
  Description:  Print regression results to dev
  Parameters:   flag: Summary: 1=verbal-header  2=math-header  4=data
                      Details: 8=verbal-header 16=math-header 32=data
  Remark:       y,B,J,M1,S,V =0 allowed if flag<8 (i.e. only summary)
 *-MH---------------------------------------------------------------------------*/
void prPCRResults(FILE *dev, double y[], int n, double lC[], int kmax, int kml, double lE, double ll, double ElogP, double VlogP,
                  double B[], int J[], double M1[], double Mt[], double S[], double V[], int flag)
/*------------------------------------------------------------------------------*/
{
  int t;

  if (flag&1) fprintf(dev,"\nFunction     \tColumn    \tKey      \t1st-row  \tDataSize \tMethod (prior) \tGlobal-Mean \tGlobal-StdDev \tSegment-StdDev \tEvidence  \tLikelihood     \tExp.Likehood\tStd.Likehood    \tRel.Likehood      \tOpt#Segments \tConfidence \t           \t           \t     \tzero \ttwo \tComments     ");
  if (flag&2) fprintf(dev,"\ndName[dmode] \tdatacol   \tdatakey  \t         \tn        \tmore by hand   \tmglob       \tsqrt(varglob) \tsqrt(varseg)   \tlog(p(y)) \tlog(p(y|mu,t,k)\tE[log p]    \tsqrt(Var[log p])\t(ll-Exp)/sqrt(Var)\tkml          \tp(kml|y)   \tp(kml-1|y) \tp(kml+1|y) \tkmax \t0    \t2   \t             ");
  if (flag&4) fprintf(dev,"\n%s           \t%d        \t%d       \t%d       \t%d       \t%s             \t%3.12lf     \t%3.12lf        \t%3.12lf       \t%3.12lf   \t%3.12lf        \t%3.12lf     \t%3.12lf         \t%3.12lf           \t%d           \t%3.12lf    \t%3.12lf    \t%3.12lf    \t%d   \t0    \t2   \t<add by hand>",
                          dName[dmode]    , datacol   , datakey  , gcnt,     n,        rName[rmode]     , mglob       , sqrt(varglob) , sqrt(varseg)   , lE        , ll          , ElogP, sqrt(VlogP), VlogP?(ll-ElogP)/sqrt(VlogP):0, kml, lC?exp(lC[kml]):0, lC?exp(lC[max(kml-1,0)]):0, lC?exp(lC[min(kml+1,kmax)]):0, kmax);
  if ((flag&56)&&(flag&7)) fprintf(dev,"\n");
  // print data vertical, since e.g. Excel takes only n<256 columns
  if (flag&8) fprintf(dev,"\n          \t             \tData        \tBreak-Probability\tPCRegression\tBayesMean   \tBayesStdDev \tMean-StdDev \tMean+StdDev \tk           \tlog-Dim.Prob.\tzero\tBreakPoints");
  if (flag&16) fprintf(dev,"\nt=1..n   \tt+0.5        \ty[t]        \tP(break|kml)     \tE(mu_k|y)   \tE[mu_t|y)   \tVar(mu_t|y) \tE-Var       \tE+Var       \tk           \tlog(p(k|y))  \t0   \tt_k        ");
  if (flag&32)
  {
    for(t=1;t<=min(n,maxnprint);t++) 
    {
      fprintf(dev,"\n%d\t%lf\t%3.12lf\t%3.12lf\t%3.12lf\t%3.12lf\t%3.12lf\t%3.12lf\t%3.12lf\t",
        gcnt+t-1,gcnt+t-0.5,y[t],B[t],Mt[t],S[t],sqrt(V[t]),S[t]-sqrt(V[t]),S[t]+sqrt(V[t]));
      if (t<=kmax) fprintf(dev,"%d\t%3.12lf\t",t,lC[t]); else fprintf(dev,"\t");
      if (t<=kml)   fprintf(dev,"0\t%lf\t",J[t]+gcnt-0.5);  else fprintf(dev,"\t\t");
    } 
    fprintf(dev,"\n");
  }
} // End of prPCRResults()

/*------------------------------------------------------------------------------*
  Description:  Perform regression and print results to dev and dev2 
                according to flag and flag2
  Parameters:   flag(2): Summary: 1=verbal-header  2=math-header  4=data
                         Details: 8=verbal-header 16=math-header 32=data
 *-MH---------------------------------------------------------------------------*/
int PCRegPrint(int flag, int flag2)
/*------------------------------------------------------------------------------*/
{
  double *y=CrVector(nmax+1,double);
  printf("Creating or Reading Data ... ");
  // Create artificial data
  int n=CreatePCData(y,nmax,dmode); 
  printf("of size n=%d ...\n",n);
  if (n==0) { printf("No suitable data\n"); return 0; }
  // Optionally estimate global parameters from data (be careful)
  EstGlobParam(y,n); 
  double **lA0=CrMatrix(n+1,n+1,double);
  double **A10=CrMatrix(n+1,n+1,double);
  double **A20=CrMatrix(n+1,n+1,double);
  CrSegDist1(y,n,lA0,A10,A20);
  double **lL=CrMatrix(kmax+1,n+1,double);
  double **lR=CrMatrix(kmax+1,n+1,double);
  int kmx=PCEvidence(lA0,n,lL,lR,kmax); // can make kmax smaller if appropriate
  double *lC=CrVector(kmax+1,double);
  int kml;
  double lE;
  double *B=CrVector(n+1,double);
  int    *J=CrVector(kmx+1,int);
  double *S=CrVector(n+1,double);
  double *V=CrVector(n+1,double);
  PCReg(lA0,A10,A20,lL,lR,n,kmx,lC,lE,kml,B,J,S,V);
  double *M1=CrVector(kml+1,double);
  double *Mt=CrVector(n+1,double);
//  for(int k=1;k<=kml;k++) M1[k]=MeanG(y,n,J[k-1]+1,J[k],rmode);
  for(int k=1;k<=kml;k++) 
  {
    M1[k]=A10[J[k-1]][J[k]];
    for(int t=J[k-1]+1;t<=J[k];t++) Mt[t]=M1[k];
  }
  double ll,ElogP,VlogP;
  PCLikelihood(y,n,Mt,varseg,ll,ElogP,VlogP);
  printf("Print Results ...\n");
  prPCRResults(dev,y,n,lC,kmx,kml,lE,ll,ElogP,VlogP,B,J,M1,Mt,S,V,flag);
  prPCRResults(dev2,y,n,lC,kmx,kml,lE,ll,ElogP,VlogP,B,J,M1,Mt,S,V,flag2);
  free(Mt); free(M1); free(V); free(S); free(J); free(B); free(lC); 
  free(lR); free(lL); free(A20); free(A10); free(lA0); free(y);
  return n;
} // End of PCRegPrint()

/********************************************************************************/
/*                            A p p l i c a t i o n s                           */
/********************************************************************************/

/*------------------------------------------------------------------------------*
  Description:  Regression for 23 Chromosomes of 31 Patients
  Parameters:   In:  File "../../../Temp/Regression/mm-copynumber-normalized.dat"
                Out: Files "../../../Temp/Regression/PatCol[03..33].dat"
 *-MH---------------------------------------------------------------------------*/
void PCRegGenomes(void)
/*------------------------------------------------------------------------------*/
{
  dmode=4;
  rmode=1;
  char dName[99];
  
  nmax=1000;
  fprintf(dev,"\n\nRegression Summary for 23 Chromosomes of 31 Patients");
  fprintf(dev,  "\n----------------------------------------------------");
  prPCRResults(dev,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1|2);
  fflush(dev);
  for(datacol=3;datacol<=33;datacol++)
  {
    sprintf(dName,"../../../Temp/Regression/PatCol%02d.dat",datacol);
    dev2=fopen(dName,"wt");
    fprintf(dev2,"\nRegression Data of PatColumn %d",datacol);
    fprintf(dev2,"\n-----------------------------");
    prPCRResults(dev2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,8|16);
    for(datakey=1;datakey<=23;datakey++) 
    {
      printf("\nAnalysing genome %d of PatColumn %d",datakey,datacol);
      printf("\n---------------------------------\n");
//      fprintf(dev,"\nGenome %d of PatColumn %d",datakey,datacol);
      PCRegPrint(4,32);
      fflush(dev); fflush(dev2);
    }
    fprintf(dev,"\n");
    fclose(dev2);
  }
} // End of PCRegGenomes()

/*------------------------------------------------------------------------------*
  Description:  Estimate and print distribution of data
  Parameters:   (e.g. to check whether Gaussian)
 *-MH---------------------------------------------------------------------------*/
void EstPrintDensity()
/*------------------------------------------------------------------------------*/
{
  int    n=10000;
  double y[11000];
  double istart=-1;
  double irange=3;
  const int inum=50;
  double istep=irange/inum;
  int    D[inum];
  double sum=0,sum2=0;

//  n=ReadDataCol("mm-copynumber-normalized.dat",y,n,2,3);
  n=ReadDataCol("data.dat",y,n,2,3);
  for(int i=0;i<inum;i++) D[i]=0;
  for(int t=1;t<=n;t++) 
  { i=(int)((log(y[t])-istart)/istep); 
    if (i>=0&&i<inum) D[i]++; 
    sum+=log(y[t]); sum2+=log(y[t])*log(y[t]);
  }
  fprintf(dev,"\nIndex\tlog(rat)\t#/n\tGauss\tCauchy\n");  
  for(i=0;i<inum;i++) fprintf(dev,"%d\t%lf\t%lf\t%lf\t%lf\n",i,istart+i*istep,D[i]/(istep*n),
    GaussDistr(istart+i*istep,sum/n,sum2/n-sum/n*sum/n),
    CauchyDistr(istart+i*istep,sum/n,sum2/n-sum/n*sum/n));
  fflush(dev);
} // End of EstPrintDensity()

/*------------------------------------------------------------------------------*
  Description:  Plot Evidence and kml as function of segment variance varseg
  Parameters:   from varseg/range..varseg*range with ngrid gridpoints
 *-MH---------------------------------------------------------------------------*/
void EvdFctSeg(FILE *dev, double range, int ngrid)
/*------------------------------------------------------------------------------*/
{
  int i=1,n,kml;
  double lE;
  double *y=CrVector(nmax+1,double);
  n=CreatePCData(y,nmax,dmode);
  fprintf(dev,"\nEvidence and ML # Segments as Function of Segment Variance");
  EstGlobParam(y,n); 
  PCEvdKml(y,n,lE,kml);
  prPCRResults(dev,y,n,0,0,kml,0,0,0,0,0,0,0,0,0,0,1|2|4);
  fprintf(dev,"\n(Some quantities refer to estimated varglob only or are not available=0)");
  fprintf(dev,"\n\n\tsegment std\tlog(evidence))\tML#Segments");
  fprintf(dev,"\n\tsqrt(varseg)\tlog(P(y|varseg))\tkml");
  fprintf(dev,"\nEstimate\t%3.12lf\t%3.12lf\t%d",sqrt(varseg),lE,kml);
  fflush(dev);
  double vsmin=varseg/range; 
  double vsmax=varseg*range;
  for(varseg=vsmin; varseg<vsmax; varseg*=pow(range,2.0/ngrid))
  {
    printf("sqrt(varseg)=%lf\n",sqrt(varseg));
    PCEvdKml(y,n,lE,kml);
    fprintf(dev,"\n%d\t%3.12lf\t%3.12lf\t%d",i++,sqrt(varseg),lE,kml);
  }
  fprintf(dev,"\n");
  free(y); fflush(dev);
} // End of EvdFctSeg()

/*------------------------------------------------------------------------------*
  Description:  Plot Evidence and kml as function of segment variance varglob
  Parameters:   from varglob/range..varglob*range with ngrid gridpoints
 *-MH---------------------------------------------------------------------------*/
void EvdFctGlob(FILE *dev, double range, int ngrid)
/*------------------------------------------------------------------------------*/
{
  int i=1,n,kml;
  double lE;
  double *y=CrVector(nmax+1,double);
  n=CreatePCData(y,nmax,dmode);
  fprintf(dev,"\nEvidence and ML # Segments as Function of Global Variance");
  EstGlobParam(y,n); 
  PCEvdKml(y,n,lE,kml);
  prPCRResults(dev,y,n,0,0,kml,0,0,0,0,0,0,0,0,0,0,1|2|4);
  fprintf(dev,"\n(Some quantities refer to estimated varglob only or are not available=0)");
  fprintf(dev,"\n\n\tglobal std\tlog(evidence))\tML#Segments");
  fprintf(dev,"\n\tsqrt(varglob)\tlog(P(y|varglob))\tkml");
  fprintf(dev,"\nEstimate\t%3.12lf\t%3.12lf\t%d",sqrt(varglob),lE,kml);
  fflush(dev);
  double vgmin=varglob/range; 
  double vgmax=varglob*range;
  for(varglob=vgmin; varglob<vgmax; varglob*=pow(range,2.0/ngrid))
  {
    printf("sqrt(varglob)=%lf\n",sqrt(varglob));
    PCEvdKml(y,n,lE,kml);
    fprintf(dev,"\n%d\t%3.12lf\t%3.12lf\t%d",i++,sqrt(varglob),lE,kml);
  }
  fprintf(dev,"\n");
  free(y); fflush(dev);
} // End of EvdFctGlob()

/*------------------------------------------------------------------------------*
  Description:  Batch Regression for all Examples in Paper
 *-MH---------------------------------------------------------------------------*/
void Exp4Paper(void)
/*------------------------------------------------------------------------------*/
{
  dev2=fopen("../../../Temp/Regression/PCRPaper.dat","wt");

//#if 0 //already done
  fprintf(dev2,"\n=================================================");
  fprintf(dev2,"\n| Regression Results for all Examples for Paper |");
  fprintf(dev2,"\n=================================================\n");
  fflush(dev2);

  // Gauss3[_~--] (low,medium,high noise)
  srand(12345);
  nmax=100; kmax=100; dmode=5; rmode=1; 
  mglob=0; varglob=1; varseg=0.01; 
  fprintf(dev2,"\nThree Gaussian Segments [_~--] with Levels +1/-1/0 and Noise %lf",sqrt(varseg));
  PCRegPrint(0,1|2|4|8|16|32|64);
  kmax=100; mglob=0; varglob=1; varseg=0.01; // redefine, since overwritten
  srand(12345);
  EvdFctSeg(dev2,10,20);

  srand(12345);
  nmax=100; kmax=100; dmode=5; rmode=1; 
  mglob=0; varglob=1; varseg=0.1; 
  fprintf(dev2,"\nThree Gaussian Segments [_~--] with Levels +1/-1/0 and Noise %lf",sqrt(varseg));
  PCRegPrint(0,1|2|4|8|16|32|64);
  kmax=100; mglob=0; varglob=1; varseg=0.1;
  srand(12345);
  EvdFctSeg(dev2,10,20);

  srand(12345);
  nmax=100; kmax=100; dmode=5; rmode=1; 
  mglob=0; varglob=1; varseg=1; 
  fprintf(dev2,"\nThree Gaussian Segments [_~--] with Levels +1/-1/0 and Noise %lf",sqrt(varseg));
  PCRegPrint(0,1|2|4|8|16|32|64);
  kmax=100; mglob=0; varglob=1; varseg=1; 
  srand(12345);
  EvdFctSeg(dev2,10,20);

  // Cauchy3[_~--] (low,medium,high noise)
  srand(12345);
  nmax=100; kmax=100; dmode=6; rmode=2; 
  mglob=0; varglob=1; varseg=0.01; 
  fprintf(dev2,"\nThree Cauchy Segments [_~--] with Levels +1/-1/0 and Noise %lf",sqrt(varseg));
  PCRegPrint(0,1|2|4|8|16|32|64);
  kmax=100; mglob=0; varglob=1; varseg=0.01; // redefine, since overwritten
  srand(12345);
  EvdFctSeg(dev2,10,20);

  srand(12345);
  nmax=100; kmax=100; dmode=6; rmode=2; 
  mglob=0; varglob=1; varseg=0.1; 
  fprintf(dev2,"\nThree Cauchy Segments [_~--] with Levels +1/-1/0 and Noise %lf",sqrt(varseg));
  PCRegPrint(0,1|2|4|8|16|32|64);
  kmax=100; mglob=0; varglob=1; varseg=0.1;
  srand(12345);
  EvdFctSeg(dev2,10,20);

  srand(12345);
  nmax=100; kmax=100; dmode=6; rmode=2; 
  mglob=0; varglob=1; varseg=1; 
  fprintf(dev2,"\nThree Cauchy Segments [_~--] with Levels +1/-1/0 and Noise %lf",sqrt(varseg));
  PCRegPrint(0,1|2|4|8|16|32|64);
  kmax=100; mglob=0; varglob=1; varseg=1; 
  srand(12345);
  EvdFctSeg(dev2,10,20);

  // Cauchy3[_~--] with medium noise, analysed by Gauss
  srand(12345);
  nmax=100; kmax=100; dmode=6; rmode=1; 
  mglob=0; varglob=1; varseg=0.1; 
  fprintf(dev2,"\nThree Cauchy Segments [_~--] with Levels +1/-1/0 and Noise %lf but analyzed with Gauss",sqrt(varseg));
  PCRegPrint(0,1|2|4|8|16|32|64);
  kmax=100; mglob=0; varglob=1; varseg=0.1;
  srand(12345);
  EvdFctSeg(dev2,10,20);

  // Gauss3[_~--] with medium noise, analysed by Cauchy
  srand(12345);
  nmax=100; kmax=100; dmode=5; rmode=2; 
  mglob=0; varglob=1; varseg=0.1; 
  fprintf(dev2,"\nThree Cauchy Segments [_~--] with Levels +1/-1/0 and Noise %lf but analyzed with Gauss",sqrt(varseg));
  PCRegPrint(0,1|2|4|8|16|32|64);
  kmax=100; mglob=0; varglob=1; varseg=0.1;
  srand(12345);
  EvdFctSeg(dev2,10,20);
//#endif

  // Patient1=Col3 Chromosome 1
  nmax=1000; kmax=100; dmode=4; datacol=3; datakey=1; rmode=1; 
  fprintf(dev2,"\nPatient%d=Col%d Chromosome %d",datacol-2,datacol,datakey);
  PCRegPrint(0,1|2|4|8|16|32|64);
  EvdFctSeg(dev2,10,20);

  // Patient3=Col5 Chromosome 9
  nmax=1000; kmax=100; dmode=4; datacol=5; datakey=9; rmode=1; 
  fprintf(dev2,"\nPatient%d=Col%d Chromosome %d",datacol-2,datacol,datakey);
  PCRegPrint(0,1|2|4|8|16|32|64);
  EvdFctSeg(dev2,10,20);

  fprintf(dev2,"\n");
  fclose(dev2);
} // End of Exp4Paper()

/*------------------------------------------------------------------------------*
  Description:  Test
 *-MH---------------------------------------------------------------------------*/
void Test(void)
/*------------------------------------------------------------------------------*/
{
  int i=1,n,kml,nmax=100;
  double lE;
  double *y=CrVector(nmax+1,double);
  srand(12345);
  printf("Test\n");
  printf("\nTest");
  fprintf(dev,"\n\tsegment std\tlog(evidence))\tML#Segments");
  fprintf(dev,"\n\tsqrt(varseg)\tlog(P(y|varseg))\tkml");
  
  for(i=1;i<30;i++)
  {
    dmode=5; rmode=1; 
    mglob=0; varglob=1; varseg=1; 
    n=CreatePCData(y,nmax,dmode);
    EstGlobParam(y,n); 
    PCEvdKml(y,n,lE,kml);
    fprintf(dev,"\n%d\t%3.12lf\t%3.12lf\t%d",i,sqrt(varseg),lE,kml);
    printf("\n%d\t%3.12lf\t%3.12lf\t%d",i,sqrt(varseg),lE,kml);
    double lEbak=lE;
    rmode=2;
    EstGlobParam(y,n); 
    PCEvdKml(y,n,lE,kml);
    fprintf(dev,"\t%3.12lf\t%3.12lf\t%d\t%lf",sqrt(varseg),lE,kml,lE-lEbak);
    printf("\t%3.12lf\t%3.12lf\t%d\t%lf",sqrt(varseg),lE,kml,lE-lEbak);
    fflush(dev);
  }
  fprintf(dev,"\n");
  free(y); fflush(dev);
} // End of Test()

/*------------------------------------------------------------------------------*
  Description:  Test2
 *-MH---------------------------------------------------------------------------*/
void Test2(void)
/*------------------------------------------------------------------------------*/
{
  nmax=1000;
  int kml;
  double lE;
  double *y=CrVector(nmax+1,double);
  srand(12345);
  printf("Test\n");
  printf("\nTest");
  fprintf(dev,"\n\tsegment std\tlog(evidence))\tML#Segments");
  fprintf(dev,"\n\tsqrt(varseg)\tlog(P(y|varseg))\tkml");
  nmax=1000; dmode=4; datacol=3; datakey=1; rmode=2; 
  int n=CreatePCData(y,nmax,dmode);
  EstGlobParam(y,n); 
  PCEvdKml(y,n,lE,kml);
  fprintf(dev,"\n%d\t%3.12lf\t%3.12lf\t%d",0,sqrt(varseg),lE,kml);
  printf("\n%d\t%3.12lf\t%3.12lf\t%d",0,sqrt(varseg),lE,kml);
  fprintf(dev,"\n");
  free(y); fflush(dev);
} // End of Test2()

/********************************************************************************/
/*                          M a i n   F u n c t i o n                           */
/********************************************************************************/

/*------------------------------------------------------------------------------*
  Description: main()
  Remark:      Activate the function you want to call
 *-MH---------------------------------------------------------------------------*/
int main(int argc, char* argv[])
/*------------------------------------------------------------------------------*/
{
  srand(12345);
  printf("PCReg started!\n");
  dev=fopen("../../../Temp/Regression/PCReg.dat","at");   // "at"=append, "wt"=overwrite
  if (!dev) panicf("Can't open PCReg.dat\n");
//  Exp4Paper(); exit(0);
//  Test();
//  PCRegGenomes(); exit(0);
//  EvdFctSeg(10,20); exit(0);
//  EvdFctGlob(100,20); exit(0);
  nmax=100; kmax=100; dmode=5; rmode=1; mglob=0; varglob=1; varseg=0.1; 
//  nmax=1000; kmax=100; dmode=4; datacol=3; datakey=1; rmode=1; 
  PCRegPrint(1|2|4|8|16|32|64,0); exit(0);
  fclose(dev);

  return 0;
}

/*------------------------------End-of-pcreg.cpp--------------------------------*/