Revision e034dc0d6fe8d41a828796f07791d50011b6bb04 authored by MoorjaniLab on 09 May 2022, 23:55:41 UTC, committed by GitHub on 09 May 2022, 23:55:41 UTC
1 parent b8f902d
fitexp2.c
#include <nag.h>
#include <nag_stdlib.h>
#include <stdio.h>
#include <limits.h>
#include <math.h>
#include <nicklib.h>
#include <nage04.h>
#include "regsubs.h"
extern int verbose ;
double *zans, *zfit, *zexp ;
static double *xxval = NULL ;
static int mval, nval, affmode = YES ;
double scorit(double *x, double *vv, int m, int n, int mode) ;
double p2s (double p) ;
double s2p (double s) ;
void loadval(double *x, int m, int n, int mode) ;
void nagopt(double *params, int n) ;
extern long seed ;
double fitexp2(double *xval, double *xfit, double *xexp, double *xco,
int m, int n, int mode, int inititer)
{
int i, k, numcols, x, iter, mm, outiter = 10, inniter ;
double y, ytry, val, vbest, ybest, yinc ;
char *fname = "qq" ;
double **xx, **fit, *try, *tbest, *ww, *wbase ;
double lo = 0.0, mul = 1.0 ;
if (seed == 0) {
seed = seednum() ;
printf("seed: %ld ", seed) ; printnl() ;
SRAND(seed) ;
}
mm = 2*m ;
if (affmode) ++mm ;
ZALLOC(try, m, double) ;
ZALLOC(tbest, m, double) ;
ZALLOC(zexp, m, double) ;
ZALLOC(zans, 2*m+1, double) ;
ZALLOC(zfit, n, double) ;
ZALLOC(ww, m, double) ;
ZALLOC(wbase, m, double) ;
vbest = 1.0e6 ;
loadval(xval, m, n, mode) ;
inniter = inititer / outiter ;
inniter += outiter ;
for (iter = 1; iter <= inititer; ++iter) {
for (k=0; k<m; ++k) {
try[k] = lo + (1.0-lo) *DRAND() ;
}
vst(ww, tbest, 1.0-mul, m) ;
vst(try, try, mul, m) ;
vvp(try, try, ww, m) ;
sortit(try, NULL, m) ;
val = scorit(xval, try, m, n, affmode) ;
if (val < vbest) {
vbest = val ;
copyarr(try, tbest, m) ;
printf("%12.6f \n", vbest*1.0e6) ;
printmat(try, 1, m) ;
}
x = iter % inniter ;
if (x == 0) mul *= 0.5 ;
}
printmat(tbest, 1, m) ;
verbose = NO ;
nagopt(tbest, m) ;
verbose = NO ;
val = scorit(xval, tbest, m, n, affmode) ;
printf("%12.6f \n", val*1.0e6) ;
printmatl(tbest, 1, m) ;
copyarr(zexp, xexp, m) ;
copyarr(zans, xco, mm) ;
copyarr(zfit, xfit, n) ;
free(try) ;
free(tbest) ;
free(zexp) ;
free(zans) ;
free(zfit) ;
free(ww) ;
free(wbase) ;
val += 1.0e-12 ;
return sqrt(val) ;
}
double scorit(double *x, double *vv, int m, int n, int mode)
{
int i, mm ;
double *eq, *rhs, *v2 ;
double z, y1, y2 ;
double *tt ;
mm = 2*m ;
if (mode) ++mm ;
ZALLOC(eq, mm*n, double) ;
ZALLOC(rhs, n, double) ;
ZALLOC(tt, mm, double) ;
ZALLOC(v2, mm, double) ;
copyarr(vv, v2, m) ;
vvt(v2+m, vv, vv, m) ;
vclear(tt, 1.0, mm) ;
for (i=0; i<n; i++) {
vvt(tt, tt, v2, 2*m) ;
copyarr(tt, eq+i*mm, mm) ;
rhs[i] = x[i] ;
}
regressit(tt, eq, rhs, n, mm) ;
copyarr(tt, zans, mm) ;
copyarr(vv, zexp, m) ;
y2 = 0.0 ;
for (i=0; i<n; i++) {
zfit[i] = vdot(tt, eq+i*mm, mm) ;
z = rhs[i] - zfit[i] ;
y2 += z*z ;
}
y2 /= (double) n ;
free(eq) ;
free(rhs) ;
free(tt) ;
free(v2) ;
if (verbose) {
printf("zz %15.9f ", y2) ;
printmat(vv, 1, m) ;
}
return y2 ;
}
double p2s (double p)
{
if (p==0.0) return -1.0e20 ;
if (p==1.0) return 1.0e20 ;
return log(p/(1-p)) ;
}
double s2p (double s)
{
double r ;
r = exp(s) ;
return r/(1+r) ;
}
void loadval(double *x, int m, int n, int mode)
{
if (xxval != NULL) free(xxval) ;
ZALLOC(xxval, n, double) ;
mval = m ;
nval = n ;
affmode = mode ;
copyarr(x, xxval, n) ;
}
void myfun(Integer n, double *xc, double *fc, Nag_Comm *comm)
{
int nparams = n ;
double *params ;
double y, yval, yfix ;
int i ;
ZALLOC(params, mval, double) ;
for (i=0; i<mval; i++) {
params[i] = s2p(xc[i]) ;
}
yval = 1.0e6*scorit(xxval, params, mval, nval, affmode) ;
free(params) ;
*fc = yval ;
}
void nagopt(double *params, int n)
{
double val ;
Nag_E04_Opt options;
static NagError fail;
double *ycoeffs ;
static int ncall = 0 ;
int nparams = n, i ;
double yval, oldyval, tmp ;
++ncall ;
e04xxc(&options);
options.print_level = Nag_NoPrint ;
if (verbose) options.print_level = Nag_Soln ;
options.optim_tol = 1.0e-10 ;
options.max_iter = n*2000 ;
options.list = FALSE ;
fail.print = TRUE;
ZALLOC(ycoeffs, nparams, double) ;
for (i=0; i<mval; i++) {
ycoeffs[i] = p2s(params[i]) ;
}
oldyval = yval = scorit(xxval, params, mval, nval, affmode) ;
/**
printf("nagopt ncall: %d value: %9.3f\n", ncall, yval) ;
for (i=0; i<n; i++) {
printf("%3d %12.6fn", i, params[i]) ;
}
*/
e04ccc(n, myfun, ycoeffs, &yval, &options, NAGCOMM_NULL, &fail);
for (i=0; i<mval; i++) {
params[i] = s2p(ycoeffs[i]) ;
}
yval = scorit(xxval, params, mval, nval, affmode) ;
printf("nagopt ncall: %d oldvalue: %15.9f newvalue: %15.9f\n", ncall, oldyval*1.0e6, yval*1.0e6) ;
free(ycoeffs) ;
}
Computing file changes ...