Commit 76752f3e authored by Aytekin Gel's avatar Aytekin Gel
Browse files

Update the PSUADE surrogate model and associated files

parent 840ca4bc

Case1_ParticleSettling/DeterministicCalibration/PSUADE/1_ResponseSurface/GPM/PSbuildRS_2.in

0 → 100644
+33 −0
Original line number Diff line number Diff line 
load ../../psData_OLH_n120_i6_o1_y2 # load data file

ana_expert      # psuade> analysis expert mode on

printlevel 3    # psuade> printlevel set from 0 to 3

rscheck         # Check for best data-fitted surrogate model

10               # Select Gaussian Process Model (GPM) Tong implementation as response surface method

1               # Quantity of interest no

y               # Perform cross-cross_validation

120             # Enter the number of groups to validate : (2 - 120) 120 fold cross-validation

n               # Random selection of leave-out groups ? (y or n)

sys mv RSFA_CV_err.m R2_GPM10_120F_RSFA_CV_err.m # Rename the CV error file before next one

sys mv RSFA_training_err.m R2_GPM10_RSFA_training_err.m # Rename the training error file before next one

rs2

256             # Grid resolution ? (32 - 256)

2               # Available response surface tools:  select Quadratic

6               # Enter the input for x axis (1 - 7) :

5               # Enter the input for y axis (1 - 7), not 1 : 2

y               # Set other nominal values automatically ? (y or n)

1               # Enter output number (1 - 1) :

n               #  Ymin and Ymax found = 2.037723e-01 3.421112e-01. Set lower threshold ? (y or n) :

n               # Set upper threshold ? (y or n) : n

sys mv matlabrs2.m R2_GPM10_rs2_x6x5.m    # Rename the response surface constructed

rs2

256             # Grid resolution ? (32 - 256)

2               # Available response surface tools:  select LINEAR

6               # Enter the input for x axis (1 - 7) :

3               # Enter the input for y axis (1 - 7), not 1 : 2

y               # Set other nominal values automatically ? (y or n)

1               # Enter output number (1 - 3) :

n               #  Ymin and Ymax found = 2.037723e-01 3.421112e-01. Set lower threshold ? (y or n) :

n               # Set upper threshold ? (y or n) : n

sys mv matlabrs2.m R2_GPM10_rs2_x6x3.m    # Rename the response surface constructed

sys ls -alrt R2_GPM*  # Check most one created after mv

q

Case1_ParticleSettling/DeterministicCalibration/PSUADE/1_ResponseSurface/GPM/R2_GPM10_120F_RSFA_CV_err.m

0 → 100644
+406 −0

File added.

Preview size limit exceeded, changes collapsed.

Case1_ParticleSettling/DeterministicCalibration/PSUADE/2_SensitivityAnalysis/matlabrssoboltsib.m

+5 −4
Original line number Diff line number Diff line
@@ -57,10 +57,11 @@ grid on 
box on

axis([0  nn+1 ymin ymax])

set(gca,'XTickLabel',[]);

th=text(1:nn, repmat(ymin-0.05*(ymax-ymin),nn,1),Str,'HorizontalAlignment','center','rotation',330);

set(th, 'fontsize', 12)

set(th, 'fontweight', 'bold')

th=text(1:nn, repmat(ymin-0.05*(ymax-ymin),nn,1),Str,'HorizontalAlignment','center');

set(th,'rotation',30)

set(th, 'fontsize', 10)

%set(th, 'fontweight', 'bold')

title('Total Order Sobol Indices (with bootstrap)','FontWeight','bold','FontSize',12)

ylabel('Total Order Sobol Index (Normalized)','FontWeight','bold','FontSize',12)

hold off

print -dpng R2_rssoboltsib
 
 No newline at end of file
 
print -dpdf R2_rssoboltsib
 
 No newline at end of file

Case1_ParticleSettling/DeterministicCalibration/PSUADE/3_Calibration/Deterministic/evalRSM_GPM10.c

0 → 100644
+201 −0
Original line number Diff line number Diff line 
/* ***********************************************/

/* GP interpolator from PSUADE.                  */

/* To estimate prediction uncertainty uncomment, */

/* dgetrs and the corresponding code segment.    */

/* ==============================================*/

#include <math.h>

#include <stdlib.h>

#include <stdio.h>

/*AIKE add string.h */

#include <string.h>



/*void dgetrs_(char*,int*,int*,double*,int*,int*,

               double*,int*,int*);*/

int initialize();

int finalize();

int interpolate(int,double*,double*,double*);

/*AIKE added int main to main */

/* main(int argc, char **argv) */

int main(int argc, char **argv)

{

  int    ii, iOne=1, nInps;  /*AIKE replaced i with ii */

  double X[6], Y, Std, finalY;  /*AIKE added finalY declaration */

  /*AIKE Add experimental dataset */

  double Exp[21][4] = {

       {1,  0.25, 0.088980934,  0.001},

       {2,  0.24, 0.088593673,  0.001},

       {3,  0.23, 0.088014071,  0.001},

       {4,  0.22, 0.087231143,  0.001},

       {5,  0.21, 0.0862336,    0.001},

       {6,  0.2,  0.085009842,  0.001},

       {7,  0.19, 0.08354796,   0.001},

       {8,  0.18, 0.081835729,  0.001},

       {9,  0.17, 0.079860606,  0.001},

       {10, 0.16, 0.077609728,  0.001},

       {11, 0.15, 0.075069909,  0.001},

       {12, 0.14, 0.072227635,  0.001},

       {13, 0.13, 0.069069065,  0.001},

       {14, 0.12, 0.065580023,  0.001},

       {15, 0.11, 0.061746,     0.001},

       {16, 0.1,  0.057552148,  0.001},

       {17, 0.09, 0.052983279,  0.001},

       {18, 0.08, 0.048023859,  0.001},

       {19, 0.07, 0.042658011,  0.001},

       {20, 0.06, 0.036869505,  0.001},

       {21, 0.05, 0.03064176,   0.001}

      };

  FILE   *fIn=NULL, *fOut=NULL;

  if (argc < 2)

  {

    printf("ERROR: not enough argument.\n");

    printf("ERROR: not enough argument.\n");

    exit(1);

  }

  fIn = fopen(argv[1], "r");

  if (fIn == NULL)

  {

    printf("ERROR: cannot open input file.\n");

    exit(1);

  }

  fscanf(fIn, "%d", &nInps);

  if (nInps != 5)  /*AIKE changed from 6 to 5 */

  {

    printf("ERROR - wrong nInputs.\n");

    exit(1);

  }

  /*AIKE changed from i<6 to i<5 and also changed counter from i to ii */

  for (ii=0; ii<5; ii++) fscanf(fIn, "%lg", &X[ii]);

  fclose(fIn);

  initialize();

  /*AIKE */

  /* interpolate(iOne, X, &Y, &Std);

  printf("Y = %e (stdev = %e)\n", Y, Std); */

  finalY = 0;

  for (ii = 0; ii < 21; ii++)

  {

    X[5] = Exp[ii][1];

    interpolate(iOne, X, &Y, &Std);

    finalY += (Y - Exp[ii][2]) * (Y - Exp[ii][2]);

  }

  finalize();

  if (argc >= 3)

  {

    fOut = fopen(argv[2], "w");

    if (fOut == NULL) {

      printf("ERROR: cannot open output file.\n");

      exit(1);

    }

    /*AIKE changed Y to finalY */

    fprintf(fOut," %e\n", finalY);

    fclose(fOut);

  }

}

/* ==========================================*/

/* Regression interpolation function         */

/* X[0], X[1],   .. X[m-1]   - first point

 * X[m], X[m+1], .. X[2*m-1] - second point

 * ... */

/* ==========================================*/

int    nInps, nSamples, nParams, *pivs=NULL;

double *XMeans=NULL,*XStds=NULL,YMean,YStd,*Thetas=NULL;

double *CMat=NULL,*CInvY=NULL, *XN=NULL;

int interpolate(int npts,double *X,double *Y,double *YStds)

{

  int    ss, ii, jj, kk;

  double expn, *xt, *yt, *zt, dist, ddata;

  xt = (double *) malloc(nInps*sizeof(double));

  yt = (double *) malloc(nSamples*sizeof(double));

  zt = (double *) malloc(nSamples*sizeof(double));

  for (ss = 0; ss < npts; ss++)

  {

    for (ii = 0; ii < nInps; ii++)

      xt[ii] = (X[ss*nInps+ii]-XMeans[ii])/XStds[ii];

    for (kk = 0; kk < nSamples; kk++)

    {

      expn = 0.0;

      for (ii = 0; ii < nInps; ii++)

      {

        dist = XN[kk*nInps+ii] - xt[ii];

        expn += pow(dist,2.000000e+00)/exp(2.000000e+00*Thetas[ii]);

      }

      expn *= 0.5;

      yt[kk] = exp(Thetas[nInps]) * exp(-expn);

    }

    Y[ss] = Thetas[nInps+2];

    for (kk = 0; kk < nSamples; kk++)

      Y[ss] += yt[kk] * CInvY[kk];

    Y[ss] = Y[ss] * YStd + YMean;

    YStds[ss] = 0.0;

    /* ==== if need to compute std dev. =====

    ddata = exp(Thetas[nInps])+exp(Thetas[nInps+1])+exp(Thetas[nInps+3]);

    LUSolve(yt, zt);

    for (kk = 0; kk < nSamples; kk++)

      ddata -= yt[kk] * zt[kk];

    if (ddata < 0) ddata = - ddata;

    YStds[ss] = sqrt(ddata);

    */

  }

  free(xt); free(yt); free(zt);

}

int initialize()

{

  int    ii, jj;

  double ddata;

  char   line[1001], word[1001];

  FILE *fp = fopen("psuade_rs.info", "r");

  if (fp == NULL)

  {

    printf("Data file (psuade_rs.info) not found.\n");

    exit(1);

  }

  while (1)

  {

    fgets(line, 500, fp);

    sscanf(line, "%s",word);

    if (!strcmp(word, "PSUADE_BEGIN")) break;

  }

  fscanf(fp, "%d %d", &nSamples, &nInps);

  XMeans = (double *) malloc(nInps*sizeof(double));

  XStds  = (double *) malloc(nInps*sizeof(double));

  for (ii = 0; ii < nInps; ii++)

  {

    fscanf(fp, "%lg", &ddata);

    XMeans[ii] = ddata;

    fscanf(fp, "%lg", &ddata);

    XStds[ii] = ddata;

  }

  fscanf(fp, "%lg", &YMean);

  fscanf(fp, "%lg", &YStd);

  fscanf(fp, "%d", &nParams);

  Thetas = (double *) malloc(nParams*sizeof(double));

  for (ii = 0; ii < nParams; ii++)

  {

    fscanf(fp, "%lg", &ddata);

    Thetas[ii] = ddata;

  }

  XN = (double *) malloc(nSamples*nInps*sizeof(double));

  for (jj = 0; jj < nSamples; jj++)

    for (ii = 0; ii < nInps; ii++)

      fscanf(fp, "%lg", &XN[jj*nInps+ii]);

  CInvY = (double *) malloc(nSamples*sizeof(double));

  for (ii = 0; ii < nSamples; ii++)

    fscanf(fp, "%lg", &CInvY[ii]);

  CMat=(double*) malloc(nSamples*nSamples*sizeof(double));

  for (jj = 0; jj < nSamples; jj++)

    for (ii = 0; ii < nSamples; ii++)

      fscanf(fp, "%lg ", &CMat[jj+ii*nSamples]);

  pivs = (int *) malloc(nSamples*sizeof(int));

  for (ii = 0; ii < nSamples; ii++)

    fscanf(fp, "%d", &pivs[ii]);

}

/* ==========================================*/

int finalize()

{

  if (XMeans != NULL) free(XMeans);

  if (XStds  != NULL) free(XStds);

  if (Thetas != NULL) free(Thetas);

  if (CInvY  != NULL) free(CInvY);

  if (CMat   != NULL) free(CMat);

  if (XN     != NULL) free(XN);

  if (pivs   != NULL) free(pivs);

}