commit 91f87cafa0ffa7e49572bfbd13092fefb1f9d7d2 Author: 大窪 貴洋 Date: Tue Jun 30 08:59:47 2026 +0900 Import existing OSILAP source code diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..7fc0b14 --- /dev/null +++ b/.gitignore @@ -0,0 +1,21 @@ +* + +!*/ + +!*.c +!*.h + +# Makefile +!Makefile +!Makefile.* +!GNUmakefile + +!README +!README.* +!LICENSE +!COPYING +!AUTHORS +!CHANGELOG +!NEWS + +!.gitignore diff --git a/BRDnewton_modify.c b/BRDnewton_modify.c new file mode 100644 index 0000000..cb7920d --- /dev/null +++ b/BRDnewton_modify.c @@ -0,0 +1,147 @@ +#include "ilt.h" +/* ####################################################################### */ +/* BRD法で最適化 */ +/* 行列サイズは以下のとおり */ +/* kz: M*N (const) */ +/* mz: M*1 (const) */ +/* KH: M*N kz*G */ +/* G: M*M (diagonal) */ +/* aI: M*M alpha*(diagonal) */ +/* ppchi: M*M */ +/* cr: M*1 */ +/* w1: M*1 */ +/* pchi: M*1 */ +/* H: N*N (diagonal) */ +/* h: N*1 = kz*cr */ +/* ####################################################################### */ +double BRDnewton_modify(const double *kz, const double *mz, + const double *aI, double *cr, + double *H, HEAD *head, const int *m, const int *n){ + const int incr = 1; + const int mm = (*m)*(*m); + const double zero = 0.0; + const double one = 1.0; + const double one_ = -1.0; + int i, info, step, m_step = 1; + int *ipiv, lwork = 4*(*m); + double gamma=1.0, pchi_norm, chi, chi_p=1e300; + double *p_cr, *p_pchi, *p_ppchi; + double mz_norm = dnrm2_(m, mz, &incr); + double *KH, *G, *ppchi; + double *w1, *dcr, *pchi, *h; + + ipiv= malloc(sizeof(int)*(*n)); + KH = malloc(sizeof(double)*(*m)*(*n)); + G = malloc(sizeof(double)*(*m)*(*m)); + ppchi = malloc(sizeof(double)*(*m)*(*m)); + w1 = malloc(sizeof(double)*(*m)); + dcr = malloc(sizeof(double)*(*m)); + pchi = malloc(sizeof(double)*(*m)); + h = malloc(sizeof(double)*(*n)); + p_cr = malloc(sizeof(double)*(*m)); + p_pchi = malloc(sizeof(double)*(*m)); + p_ppchi = malloc(sizeof(double)*(*m)*(*m)); + + double *work, *wi, *wr, *vl, *vr, lambda; + work = malloc(sizeof(double)*lwork); + wr = malloc(sizeof(double)*(*m)); + wi = malloc(sizeof(double)*(*m)); + vl = malloc(sizeof(double)*(*m)*(*m)); + vr = malloc(sizeof(double)*(*m)*(*m)); + + for (step=0;; step++){ + /* hi = np.dot(k0.T, cr) */ + dgemv_("T", m, n, &one, kz, m, cr, &incr, &zero, h, &incr); + /* H (diagonal)の計算 */ + for (i=0; i<*n; i++){ + if (h[i]>0) H[i*(*n+1)]=1.0; + else H[i*(*n+1)]=0.0; + } + /* Gr = kz*H*kzT */ + dgemm_("N", "N", m, n, n, &one, kz, m, H, n, &zero, KH, m, &info); + dgemm_("N", "T", m, m, n, &one, KH, m, kz, m, &zero, G, m, &info); + + /* chi=0.5*cr*(G+aI)*cr - cr*mz */ + daxpy_(&mm, &one, aI, &incr, G, &incr); /* G = G+aI */ + dcopy_(&mm, G, &incr, ppchi, &incr); /* G = chi'' */ + dgemv_("N", m, m, &one, G, m, cr, &incr, &zero, w1, &incr);/* w1=(G+aI)cr */ + chi = 0.5*ddot_(m, cr, &incr, w1, &incr); + chi = chi - ddot_(m, cr, &incr, mz, &incr); + /* chi' = (G+alphaI)cr-mz */ + dcopy_(m, mz, &incr, pchi, &incr); /* cr' = mz */ + dscal_(m, &one_, pchi, &incr); /* cr' = -mz */ + daxpy_(m, &one, w1, &incr, pchi, &incr); /* cr' = w1 -mz */ + pchi_norm = dnrm2_(m, pchi, &incr); + + /* chi'とchi''からdcr=w1=-ppchi/pchiを求める */ + if (chi_p > chi){ + dcopy_(m, cr, &incr, p_cr, &incr); /* crを保存しておく */ + dcopy_(&mm, ppchi, &incr, p_ppchi, &incr); /* chi''を保存しておく */ + dcopy_(m, pchi, &incr, p_pchi, &incr); /* chi'を保存しておく */ + dgesv_(m, &incr, ppchi, m, ipiv, pchi, m, &info); /* phiにdcrが入る */ + dcopy_(m, pchi, &incr, dcr, &incr); /* dcrを保存しておく */ + gamma = -1; + daxpy_(m, &gamma, dcr, &incr, cr, &incr); /* cr - dcr */ + } + else { + dcopy_(m, p_cr, &incr, cr, &incr); /* crを元に戻す */ + dcopy_(m, p_pchi, &incr, pchi, &incr); /* cr'を元に戻す */ + dcopy_(&mm, p_ppchi, &incr, ppchi, &incr); /* cr''を元に戻す */ + /* dcr方向へのステップを減らしてcrを作る */ + gamma = -pow(0.5, (double)m_step); + daxpy_(m, &gamma, dcr, &incr, cr, &incr); /* cr - dcr */ + m_step++; + /* 100ステップでchiが減少しなかったら越えたらppchiの正定値を疑う */ + /* Levenberg-Marquardtの方法。固有値を求めて正定値にする */ + if (m_step > 10000){ + printf("Levenberg-Marquardt...\n"); + dgeev_("N","N", m, ppchi, m, wr, wi, vl, m, vr, m, work, &lwork, &info); + dcopy_(&mm, p_ppchi, &incr, ppchi, &incr); /* cr''を元に戻す */ + lambda = dnrm2_(&mm, ppchi, &incr); + lambda = wr[idamax_(m, wr, &incr)-1] * lambda; + if (lambda < 0) lambda = -lambda; + for (i=0; i<*m; i++) ppchi[i*(*m+1)] = ppchi[i*(*m+1)] + lambda; + dgesv_(m, &incr, ppchi, m, ipiv, pchi, m, &info); /* phiにdcrが入る */ + dcopy_(m, pchi, &incr, dcr, &incr); + m_step = 1; + } + continue; + } + gamma = 1.0; + m_step = 0; + chi_p = chi; + /* 途中経過を表示 */ + printf("%4d: ", step+1); + /* printf("Dchi: %10.3e", delta_chi); */ + printf(" chi: %10.4g", chi); + printf(" chi': %10.4g", pchi_norm); + /* printf(" res : %8.3e", dnrm2_(m, w1, &incr)); */ + printf(" alpha: %10.4e", aI[0]); + printf(" loop: %d\n", head->alpha_loop+1); + fflush(stdout); + /* cr'で収束チェック */ + /* if (dnrm2_(m, w1, &incr) < head->newton_tol) break; */ + if (pchi_norm/mz_norm < head->newton_tol) { + dcopy_(m, mz, &incr, w1, &incr); + dgesv_(m, &incr, p_ppchi, m, ipiv, w1, m, &info); /* phiにdcrが入る */ + return dnrm2_(m, w1, &incr); + } + if (step == head->newton_loop_max) { + printf("Iteration reached maximum newton loop.\n"); + break; + } + } + free(work); + free(ipiv); + free(KH); + free(G); + free(ppchi); + free(w1); + free(dcr); + free(pchi); + free(h); + free(p_cr); + free(p_pchi); + free(p_ppchi); + return 0; +} diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..34d4989 --- /dev/null +++ b/Makefile @@ -0,0 +1,21 @@ +CC = gcc-mp-14 +FC = gfortran-mp-14 +CFLAGS = -O2 -Wall +LIBS = libopenblas-mac.a -lpthread \ + -static-libgfortran \ + -static-libgcc \ + -static-libquadmath +LDFLAGS = +OBJS = ilt.o nnls.o command.o output.o input.o BRDnewton_modify.o alpha_loop.o +PROGRAM = osilap + +all: $(PROGRAM) + +$(PROGRAM): $(OBJS) + $(FC) $(OBJS) $(LDFLAGS) $(LIBS) -o $(PROGRAM) + +clean: + rm -f *.o *~ $(PROGRAM) + +install: $(PROGRAM) + install -s $(PROGRAM) $(DEST) diff --git a/Makefile.icc b/Makefile.icc new file mode 100644 index 0000000..6abfb8b --- /dev/null +++ b/Makefile.icc @@ -0,0 +1,19 @@ +CC = icc +FC = icc +CFLAGS = -Wall -O3 -static-intel -qopenmp -parallel +LIBS = -mkl -lpthread -static-intel -qopenmp -qopenmp-link=static + +LDFLAGS = +OBJS = ilt.o nnls.o command.o output.o input.o BRDnewton_modify.o alpha_loop.o +PROGRAM = osilap + +all: $(PROGRAM) + +$(PROGRAM): $(OBJS) + $(FC) $(OBJS) $(LDFLAGS) $(LIBS) -o $(PROGRAM) + +clean: + rm -f *.o *~ $(PROGRAM) + +install: $(PROGRAM) + install -s $(PROGRAM) $(DEST) diff --git a/Makefile.win32 b/Makefile.win32 new file mode 100644 index 0000000..337327b --- /dev/null +++ b/Makefile.win32 @@ -0,0 +1,19 @@ +CC = i686-w64-mingw32-gcc +FC = i686-w64-mingw32-gfortran +CFLAGS = -O2 -Wall +LIBS = libopenblas-win32.a -lpthread -static-libgcc -Wl,-Bstatic -lgfortran -lquadmath -Wl,-Bdynamic -lm + +LDFLAGS = +OBJS = ilt.o nnls.o command.o output.o input.o BRDnewton_modify.o alpha_loop.o +PROGRAM = osilap + +all: $(PROGRAM) + +$(PROGRAM): $(OBJS) + $(FC) $(OBJS) $(LDFLAGS) $(LIBS) -o $(PROGRAM) + +clean: + rm -f *.o *~ $(PROGRAM) + +install: $(PROGRAM) + install -s $(PROGRAM) $(DEST) diff --git a/alpha_loop.c b/alpha_loop.c new file mode 100644 index 0000000..89b98d2 --- /dev/null +++ b/alpha_loop.c @@ -0,0 +1,254 @@ +#include "ilt.h" + +/* ####################################################################### */ +/* 最適化したfから圧縮無しのσ^2を求める */ +/* ####################################################################### */ +double CalcSigma(MATRIX *K1, MATRIX *K2, double *f, HEAD *head, INDATA *d){ + int i, incr = 1; + int m1 = K1->m, m2 = K2->m; + int n1 = K1->n, n2 = K2->n; + int nn = n1*n2, mm = m1*m2; + int info; + double one=1.0, one_=-1.0, zero = 0.0, chi, sigma=0; + double *KF, *F, *M; + div_t temp; + + if (head->dimension == 2){ + KF = malloc(sizeof(double)*m1*n2); + F = malloc(sizeof(double)*nn); + M = malloc(sizeof(double)*m1*m2); + /* fから行列Fを生成 */ + for (i=0; ia, &m1, + F, &n1, &zero, KF, &m1, &info); + /* (K1*F)*K2Tの計算 */ + dgemm_("N", "T", &m1, &m2, &n2, &one, KF, &m1, + K2->a, &m2, &zero, M, &m1, &info); + /* M = -Mexp + Mcalの計算 */ + daxpy_(&mm, &one_, d->m.a, &incr, M, &incr); + /* ||Mcal - Msim|| */ + chi = dnrm2_(&mm, M, &incr); + /* sigma = σ^2 = ||KF-M||^2/m */ + sigma = chi*chi/mm; + free(KF); + free(F); + free(M); + } + if (head->dimension == 1){ + M = malloc(sizeof(double)*m1); + dgemv_("N", &m1, &n1, &one, K1->a, &m1, f, &incr, &zero, M, &incr); + /* M = -Mexp + Mcalの計算 */ + daxpy_(&m1, &one_, d->m.a, &incr, M, &incr); + chi = dnrm2_(&m1, M, &incr); + /* sigma = σ^2 = ||KF-M||^2/m */ + sigma = chi*chi/m1; + free(M); + } + return sigma; +} + +/* ####################################################################### */ +/* Scurve法でalpha_optを求める */ +/* ####################################################################### */ +void ScurveAlphaLoop(MATRIX *K0, MATRIX *Mz, HEAD *head, double *f, + MATRIX *k1, MATRIX *k2, INDATA *d){ + int m = K0->m; + int n = K0->n; + int i, j, pnt = 0, incr = 1; + int alpha_num=head->alpha_num; + double *kz, *mz, *cr, *aI, *H; + double alpha, sigma; + double zero = 0.0, one = 1.0; + double alpha_s = head->alpha_start, alpha_e = head->alpha_end; + double s_target=head->s_taget, alpha_log; + double *alpha_list, *sigma_list, *grad; + double *t = NULL; + + kz = K0->a; + mz = Mz->a; + cr = malloc(sizeof(double)*m); + aI = malloc(sizeof(double)*m*m); + H = malloc(sizeof(double)*n*n); + alpha_list = (double *)malloc(sizeof(double)*alpha_num); + sigma_list = (double *)malloc(sizeof(double)*alpha_num); + grad = malloc(sizeof(double)*(alpha_num-1)); + /* alphaのグリッド生成 */ + for (i=0; i=0; j--){ + /* alphaI (diagonal)の初期化*/ + for (i=0; i0 だけが解 */ + for (i=0; ialpha_loop++; + } + + for (i=0; i s_target && pnt == 0) pnt = i; + } + double alpha_0, alpha_1, grad_0, grad_1; + double sigma_0, sigma_1; + int step; + alpha_0 = alpha_list[pnt-1]; + sigma_0 = sigma_list[pnt-1]; + alpha_1 = alpha_list[pnt+1]; + sigma_1 = sigma_list[pnt+1]; + for (step=1;; step++){ + pnt = alpha_num+step; + alpha = alpha_0 + 0.5*(alpha_1 - alpha_0); + if ((t = (double *)realloc(alpha_list, pnt*sizeof(double))) != NULL){ + alpha_list = t; + alpha_list[alpha_num+step-1] = alpha; + } + for (i=0; ialpha = alpha; + if ((alpha_1 - alpha_0) < head->alpha_tol) break; + printf("delta alpha = %g\n", alpha_1 - alpha_0); + printf("grad0 = %g\n", grad_0); + printf("grad1 = %g\n", grad_1); + head->alpha_loop++; + } + + head->alpha_list = malloc(pnt*sizeof(double)); + head->sigma_list = malloc(pnt*sizeof(double)); + head->alpha_num = pnt; + for (i=0; ialpha_list[i] = alpha_list[i]; + head->sigma_list[i] = sigma_list[i]; + } + + free(cr); + free(aI); + free(H); + free(alpha_list); + free(sigma_list); + free(grad); +} + +void BRDAlphaLoop(MATRIX *K0, MATRIX *Mz, HEAD *head, double *f, + MATRIX *k1, MATRIX *k2, INDATA *d){ + int m = K0->m; + int n = K0->n; + int i, incr = 1; + double *kz, *mz; + double *cr, *cr_n; + double *aI, *H; + double zero = 0.0, one = 1.0, one_ = -1.0; + double alpha_opt, cr_norm, sigma, sigma_z; + + if (head->alpha0 == 0) { + head->alpha = 0.0; + return; + } + + kz = K0->a; + mz = Mz->a; + cr = malloc(sizeof(double)*m); + cr_n = malloc(sizeof(double)*m); + aI = malloc(sizeof(double)*m*m); + H = malloc(sizeof(double)*n*n); + + /* 初期のfからcrを求める cr=-(Kf-m)/alpha */ + dgemv_("N", &m, &n, &one, kz, &m, f, &incr, &zero, cr, &incr); + daxpy_(&m, &one_, mz, &incr, cr, &incr); + alpha_opt = -1.0/head->alpha0; + dscal_(&m, &alpha_opt, cr, &incr); + head->alpha = head->alpha0; + + /* alphaI (diagonal) とH (diagonal)の初期化 */ + for (i=0; i<(n*n); i++) H[i] = 0.0; + for (i=0; ialpha_loop=0;; head->alpha_loop++){ + /* alphaI (diagonal)の初期化*/ + for (i=0; ialpha; + /* BRDでcrの最適化 */ + cr_norm = BRDnewton_modify(kz, mz, aI, cr, H, head, &m, &n); + /* crからfの計算f f= kz'*cr, fi>0 だけが解 */ + dgemv_("T", &m, &n, &one, kz, &m, cr, &incr, &zero, f, &incr); + for (i=0; isigma0 != 0.0) sigma = head->sigma0; + printf("c_norm = %g\n", cr_norm); + + switch(head->Atype){ + case 0: + alpha_opt = sigma/sigma_z; + break; + case 1: + alpha_opt = sigma_z/dnrm2_(&n, f, &incr)/(sqrt(1-sqrt(2*m)/m)); + break; + case 2: + alpha_opt = sigma_z/dnrm2_(&n, f, &incr); + break; + case 3: + alpha_opt = sigma/dnrm2_(&n, f, &incr)/(sqrt(1-sqrt(2*m)/m)); + break; + } + + /* alphaの収束チェック */ + if (fabs(alpha_opt - head->alpha)/head->alpha < head->alpha_tol) + break; + /* alphaをmixingしてupdate */ + if (head->alpha_loop > 1) + head->alpha = alpha_opt; + else + head->alpha_loop++; + if (head->alpha_loop == head->alpha_loop_max) { + printf("Iteration reached maximum alpha loop.\n"); + break; + } + } + free(cr_n); + free(cr); + free(aI); + free(H); +} diff --git a/command.c b/command.c new file mode 100644 index 0000000..4c332a3 --- /dev/null +++ b/command.c @@ -0,0 +1,195 @@ +#include +#include +#include "ilt.h" + +void Error(){ + printf("Usage: osilap [options] infile.int\n\n"); + printf("Options:\n"); + printf(" -i confirm calculation condition\n"); + printf(" -h --help show this help message\n"); + printf(" --dimension=%%d one- or two-dimension [2]\n"); + printf(" --alpha=%%f initial value of alpha (only BRD) [0.1]\n"); + printf(" --sigma=%%f standard error, 0=auto [0.0]\n"); + printf(" --sfactor=%%f sfactor for SVD [1e-4]\n"); + printf(" --output=%%s basename of output file [basename]\n"); + printf(" --TXspace=%%s TX space (log | linear) [log]\n"); + printf(" --TXtype=%%s TX data type (T1 | T1inf | T1p | T2 | D) [T1]\n"); + printf(" --TXmin=%%f TX minimum [0.0001]\n"); + printf(" --TXmax=%%f TX maximum [10]\n"); + printf(" --TXnum=%%d TX points [50]\n"); + printf(" --TYspace=%%s TY space (log | linear) [log]\n"); + printf(" --TYtype=%%s TY data type (T1 | T2 | D) [T2]\n"); + printf(" --TYmin=%%f TY minimum [0.0001]\n"); + printf(" --TYmax=%%f TY maximum [10]\n"); + printf(" --TYnum=%%d TY points [50]\n"); + printf(" --alpha_loop_max=%%d maximum iteration for alpha [5000]\n"); + printf(" --alpha_tol=%%f alpha tolerance [1e-3]\n"); + printf(" --newton_loop_max=%%d maximum iteration for newton method [5000]\n"); + printf(" --newton_tol=%%f newton tolerance [1e-8]\n"); + printf(" --Scurve=%%f target for dlog10(chi)/dlog10(alpha) [0.0]\n"); + printf(" --Smin=%%f search grid for s-curve [1e+6]\n"); + printf(" --Smax=%%f search grid for s-curve [1e-6]\n"); + printf(" --Snum=%%d search grid number for s-curve [20]\n"); + printf(" --Atype=%%d calculation method of alpha for BRD routine [0]\n"); + exit(1); +} + +void Confirm(HEAD *head, INDATA *d){ + int ch; + printf("\n"); + printf("-------- Solving condition -------\n"); + printf(" dimension: %d\n", head->dimension); + printf(" sfactor: %g\n", head->sfactor); + printf(" initial alpha: %g\n", head->alpha0); + printf(" sigma: %g\n", head->sigma0); + printf(" alpha max Loop: %d\n", head->alpha_loop_max); + printf(" alpha tol: %g\n", head->alpha_tol); + printf(" Newton max Loop: %d\n", head->newton_loop_max); + printf(" Newton tol: %g\n", head->newton_tol); + printf(" A_opt method: %d\n", head->Atype); + printf("--------- TX information ---------\n"); + printf(" Grid space: %s\n", head->TXspace); + printf(" Data type: %s\n", head->TXtype); + printf(" (TXmin, TXmax, TXpoints) ="); + printf(" (%g, %g, %d)\n", head->TXmin, head->TXmax, head->TXnum); + if (head->dimension>1){ + printf("--------- TY information ---------\n"); + printf(" Grid space: %s\n", head->TYspace); + printf(" Data type: %s\n", head->TYtype); + printf(" (TYmin, TYmax, TYpoints) ="); + printf(" (%g, %g, %d)\n", head->TYmin, head->TYmax, head->TYnum); + } + if (head->s_taget != 0){ + printf("------ S-curve configuration -----\n"); + printf(" Alpha grid space: \n"); + printf(" (start, end, points) ="); + printf(" (%g, %g, %d)\n", head->alpha_start, head->alpha_end, head->alpha_num); + } + printf("----------- Input data -----------\n"); + printf(" Input file: %s\n", head->infile); + printf(" Input data size:"); + if (head->dimension == 1){ + printf(" %d\n", d->u.m); + } + if (head->dimension == 2){ + printf(" (X, Y)=(%d, %d)\n", d->u.m, d->v.m); + } + printf(" Output basename: %s\n", head->outbase); + printf("\n"); + + fflush(stdout); + if (head->confirm == 1){ + printf("Start calculation? [y/n]: "); + while ((ch = getchar()) != EOF) { /* 文字の入力 */ + if (ch == 'y') break; + if (ch == 'n') exit(1); + } + } + fflush(stdin); +} + +/* ####################################################################### */ +/* ベース名を決める */ +/* ####################################################################### */ +void GetBaseName(HEAD *head) +{ + int nLen; + char *res = NULL; + char szPath[1024]; + + strcpy(szPath, head->infile); + nLen = sizeof (szPath); + memset(head->outbase, 0, sizeof(head->outbase)); + if ((nLen > 0) && (nLen <= 1024)) { + res = strrchr( szPath, '.' ); + if(res != NULL) + *res = '\0'; + } + strcpy(head->outbase, szPath); +} + + +int SetArgment(int argc, char * argv[], HEAD *head){ + int opt, option_index; + struct option long_options[] = { + {"help", no_argument, NULL, 'h'}, + {"confirm", no_argument, NULL, 'i'}, + {"DX", required_argument, NULL, 'A'}, + {"DY", required_argument, NULL, 'B'}, + {"TXtype", required_argument, NULL, 'C'}, + {"TXspace", required_argument, NULL, 'D'}, + {"TXmin", required_argument, NULL, 'E'}, + {"TXmax", required_argument, NULL, 'F'}, + {"TXnum", required_argument, NULL, 'G'}, + {"TYtype", required_argument, NULL, 'H'}, + {"TYspace", required_argument, NULL, 'I'}, + {"TYmin", required_argument, NULL, 'J'}, + {"TYmax", required_argument, NULL, 'K'}, + {"TYnum", required_argument, NULL, 'L'}, + {"output", required_argument, NULL, 'M'}, + {"sfactor", required_argument, NULL, 'N'}, + {"dimension", required_argument, NULL, 'O'}, + {"alpha", required_argument, NULL, 'P'}, + {"alpha_loop_max", required_argument, NULL, 'Q'}, + {"alpha_tol", required_argument, NULL, 'R'}, + {"newton_loop_max", required_argument, NULL, 'S'}, + {"newton_tol", required_argument, NULL, 'T'}, + {"sigma", required_argument, NULL, 'U'}, + {"Scurve", required_argument, NULL, 'V'}, + {"Smin", required_argument, NULL, 'W'}, + {"Smax", required_argument, NULL, 'X'}, + {"Snum", required_argument, NULL, 'Y'}, + {"Atype", required_argument, NULL, 'Z'}, + {0, 0, 0, 0}// 配列の最後はすべて0で埋める + }; + + //opterr = 0;/* エラーメッセージを非表示にする */ + if (argc == 1) Error(); + + //while((opt = getopt(argc, argv, "mes:")) != -1){ + while((opt = getopt_long(argc, argv, "ih", + long_options, &option_index)) != -1){ + switch(opt){ + case 'i': + head->confirm = 1; + break; + case 'A': head->DXnum = atoi(optarg); break; + case 'B': head->DYnum = atoi(optarg); break; + case 'C': strcpy(head->TXtype, optarg); break; + case 'D': strcpy(head->TXspace, optarg); break; + case 'E': head->TXmin = atof(optarg); break; + case 'F': head->TXmax = atof(optarg); break; + case 'G': head->TXnum = atoi(optarg); break; + case 'H': strcpy(head->TYtype, optarg); break; + case 'I': strcpy(head->TYspace, optarg); break; + case 'J': head->TYmin = atof(optarg); break; + case 'K': head->TYmax = atof(optarg); break; + case 'L': head->TYnum = atoi(optarg); break; + case 'M': strcpy(head->outbase, optarg); break; + case 'N': head->sfactor = atof(optarg); break; + case 'O': head->dimension = atoi(optarg); break; + case 'P': head->alpha0 = atof(optarg); break; + case 'Q': head->alpha_loop_max = atoi(optarg); break; + case 'R': head->alpha_tol = atoi(optarg); break; + case 'S': head->newton_loop_max = atoi(optarg); break; + case 'T': head->newton_tol = atof(optarg); break; + case 'U': head->sigma0 = atof(optarg); break; + case 'V': head->s_taget = atof(optarg); break; + case 'W': head->alpha_start = atof(optarg); break; + case 'X': head->alpha_end = atof(optarg); break; + case 'Y': head->alpha_num = atoi(optarg); break; + case 'Z': head->Atype = atoi(optarg); break; + // 解析できないオプションが見つかった場合は「?」を返す + // オプション引数が不足している場合も「?」を返す + case '?': + Error(); + } + } + strcpy(head->infile, argv[optind]); + if (strcmp(head->outbase, "")==0) GetBaseName(head); + /* if (head->dimension == 1 && head->TXnum < 100){ */ + /* head->TXnum = 200; */ + /* strcpy(head->TXtype, "T2"); */ + /* } */ + return 0; +} diff --git a/ilt.c b/ilt.c new file mode 100644 index 0000000..5d8ae7b --- /dev/null +++ b/ilt.c @@ -0,0 +1,522 @@ +#include "ilt.h" + +/* ####################################################################### */ +/* Setting Parameters */ +/* default values */ +/* ####################################################################### */ +void SetParameter(HEAD *head){ + head->newton_tol = 1e-8; + head->alpha_tol = 1e-4; + head->confirm = 0; + head->dimension = 2; + strcpy(head->outbase, ""); + head->Atype = 0; + head->sfactor = 1.0e-4; + head->alpha0 = 0.1; + head->alpha_loop_max = 5000; + head->newton_loop_max = 5000; + head->sigma0 = 0; + head->DXnum = 0; + head->DYnum = 0; + head->TXnum = 50; + head->TXmin = 1.0e-4; + head->TXmax = 10; + strcpy(head->TXspace, "log"); + strcpy(head->TXtype, "T1"); + + head->TYnum = 50; + head->TYmin = 1.0e-4; + head->TYmax = 10; + strcpy(head->TYspace, "log"); + strcpy(head->TYtype, "T2"); + + head->newton_loop = 0; + head->alpha_loop = 0; + + head->s_taget = 0.0; + head->alpha_start = -4; /* log10 space */ + head->alpha_end = 6; /* log10 space */ + head->alpha_num = 20; +} + +/* ####################################################################### */ +/* Transpose Matrix */ +/* ####################################################################### */ +void Transpose(MATRIX *A, MATRIX *B){ + int x; + int i, j, p; + div_t temp; + + B->m = A->n; + B->n = A->m; + B->a = (double*)malloc(sizeof(double)*B->m*B->n); + for (x=0; xm*A->n; x++){ + temp = div(x, A->m); + i = temp.rem; j = temp.quot; + p = i*B->m + j; + B->a[p] = A->a[x]; + } +} + +/* ####################################################################### */ +/* Kronecker (Tensor) product */ +/* ####################################################################### */ +void Kronecker(MATRIX *k1, MATRIX *k2, MATRIX *k0){ + int i, j, k, l; + int row, col; + int p; + int x, y; + div_t temp; + + k0->m = k1->m*k2->m; + k0->n = k1->n*k2->n; + k0->a = (double*)malloc(sizeof(double)*k0->m*k0->n); + + for (x=0; xm*k1->n; x++){ + temp = div(x, k1->m); + i = temp.rem; j = temp.quot; + for (y=0; ym*k2->n; y++){ + temp = div(y, k2->m); + k = temp.rem; l = temp.quot; + row = i*k2->m + k; + col = j*k2->n + l; + p = col * k0->m + row; + k0->a[p] = k1->a[x] * k2->a[y]; + } + } +} + +/* ####################################################################### */ +/* Lexicographic row major */ +/* ####################################################################### */ +void Lexicographic(MATRIX *A, MATRIX *Al){ + int i, p; + const int incr = 1; + const int mn = A->m*A->n; + div_t temp; + + Al->a = (double*)malloc(sizeof(double)*A->m*A->n); + Al->m = A->m * A->n; + Al->n = 1; + + dcopy_(&mn, A->a, &incr, Al->a, &incr); + for (i=0; in); + p = temp.rem * A->m + temp.quot; + Al->a[i] = A->a[p]; + } +} + +/* ####################################################################### */ +/* 行列の掛け算 C = A.a*B.a + beta*C.a */ +/* alpha=beta=1 */ +/* ####################################################################### */ +void MultiplyMatrix(MATRIX *A, MATRIX *B, MATRIX *C, + char aopt[], char bopt[]){ + int m, k, n; + int lda, ldb, ldc; + double alpha = 1.0; + double beta = 0.0; + int info; + + if (strcmp(aopt, "N")==0){ + m = A->m; k = A->n; + lda = A->m; ldb = B->m; ldc = A->m; + C->m = A->m; + } + if (strcmp(aopt, "T")==0){ + m = A->n; k = A->m; + lda = A->m; ldb = B->m; ldc = A->n; + C->m = A->n; + } + if (strcmp(bopt, "N")==0){ + n = B->n; ldb = B->m; + C->n = B->n; + } + if (strcmp(bopt, "T")==0){ + n = B->m; ldb = B->m; + C->n = B->m; + } + + C->a = (double*)malloc(sizeof(double)*C->m*C->n); + dgemm_(aopt, bopt, &m, &n, &k, &alpha, A->a, &lda, + B->a, &ldb, &beta, C->a, &ldc, &info); +} + +/* ####################################################################### */ +/* 配列の出力 */ +/* ####################################################################### */ +void PrintList(const double *x, const int m){ + int i; + + printf("Vector size: %d\n", m); + for (i=0; im = points; + x->n = 1; + x->a = (double*)malloc(sizeof(double)*points); + if (strcmp(type, "log")==0){ + start = log10(start); + end = log10(end); + dx = (end - start)/(double)(points-1); + for (i=0; ia[i] = pow(10, start+i*dx); + } + else if (strcmp(type, "linear")==0){ + start = start; + end = end; + dx = (end - start)/(double)(points-1); + for (i=0; ia[i] = start+i*dx; + } + else { + printf("Kernel space %s is strange.", type); + Error(); + } +} + +/* ####################################################################### */ +/* kernel生成 */ +/* ####################################################################### */ +void MakeKernel(MATRIX *t, MATRIX *x, char *type, MATRIX *k){ + int i, j, c=0; + + if ((strcmp(type, "T1inf") == 0) || + (strcmp(type, "T2inf") == 0) || + (strcmp(type, "T2ginf") == 0) ){ + k->a = (double*)malloc(sizeof(double)*t->m*(x->m+1)); + k->m = t->m; + k->n = x->m + 1; + } + else { + k->a = (double*)malloc(sizeof(double)*t->m*x->m); + k->m = t->m; + k->n = x->m; + } + if (strcmp(type, "T1") == 0 || + strcmp(type, "T1inf") == 0){ + for (i=0; im; i++){ + for (j=0; jm; j++){ + k->a[c] = 1.0-2.0*exp(-t->a[j]/x->a[i]); + c++; + } + } + } + if (strcmp(type, "T2") == 0 || + strcmp(type, "T2inf") == 0 || + strcmp(type, "T1p") == 0){ + for (i=0; im; i++){ + for (j=0; jm; j++){ + k->a[c] = exp(-t->a[j]/x->a[i]); + c++; + } + } + } + if (strcmp(type, "D") == 0 ){ + for (i=0; im; i++){ + for (j=0; jm; j++){ + k->a[c] = exp(-t->a[j]*x->a[i]); + c++; + } + } + } + if (strcmp(type, "T2g") == 0 || + strcmp(type, "T2ginf") == 0){ + for (i=0; im; i++){ + for (j=0; jm; j++){ + k->a[c] = exp(-0.5*(t->a[j]/x->a[i])*(t->a[j]/x->a[i])); + c++; + } + } + } + + /* base lineの継ぎ足し */ + if (strcmp(type, "T1inf") == 0 || + strcmp(type, "T2inf") == 0|| + strcmp(type, "T2ginf") == 0){ + for (j=0; jm; j++){ + k->a[c] = 1.0; + c++; + } + } +} + +/* ####################################################################### */ +/* SVD compress */ +/* A = U * S * Vt */ +/* ####################################################################### */ +void SvdCompress(MATRIX *A, MATRIX *uz, MATRIX *sz, MATRIX *vz, + double sfactor){ + int i, j; + const int incr = 1; + const int m = A->m, n= A->n, mn = A->m*A->n; + const int lwork=5*m*n; + int info; + double *U, *Vt, *Aa, *work; + MATRIX S; + + (m < n) ? (S.m = A->m) : (S.m = A->n); + S.n = 1; + Aa = malloc(sizeof(double)*m*n); + U = malloc(sizeof(double)*m*m); + S.a = malloc(sizeof(double)*m*n); + Vt = malloc(sizeof(double)*n*n); + work = malloc(sizeof(double)*lwork); + dcopy_(&mn, A->a, &incr, Aa, &incr); + dgesvd_("A", "A", &m, &n, Aa, &m, S.a, + U, &m, Vt, &n, work, &lwork, &info); + if (info != 0){ + printf("svd is fault!!! info=%d\n", info); + exit(1); + } + + /* 切り捨てずに採用するランク数の計算 */ + for (i=0; S.a[i]>S.a[0]*sfactor; i++); + sz->m = i; sz->n = i; + sz->a = (double*)malloc(sizeof(double) * sz->m * sz->n); + uz->m = m; uz->n = i; + uz->a = (double*)malloc(sizeof(double) * uz->m * uz->n); + vz->m = i; vz->n = n; + vz->a = (double*)malloc(sizeof(double) * vz->m * vz->n); + + /* 圧縮した行列の生成 */ + for (i=0; i < (sz->m*sz->n); i++) sz->a[i] = 0.0; /* 初期化 */ + for (i=0; i < sz->n; i++) sz->a[i*(sz->n+1)] = S.a[i]; /* diagonalだけ */ + for (i=0; i < (uz->n*uz->m); i++) uz->a[i] = U[i]; + for (i=0; i < (vz->n); i++){ + for (j=0; j < vz->m; j++){ + vz->a[j+i*vz->m] = Vt[j+i*n]; + } + } + free(S.a); + free(Aa); + free(U); + free(Vt); + free(work); +} + +/* ####################################################################### */ +/* NNLSで初期値を決める */ +/* ####################################################################### */ +void Calcf0(MATRIX *k0, MATRIX *mz, double *f0){ + double rnorm; + const int incr = 1; + const int m=k0->m; + const int n=k0->n; + const int mn=k0->m*k0->n; + int index[k0->n]; + int mode; + double *mn_mem, *m_mem, *m_mem2, *n_mem; + + mn_mem = malloc(sizeof(double)*k0->m*k0->n); + m_mem = malloc(sizeof(double)*k0->m); + m_mem2 = malloc(sizeof(double)*k0->m); + n_mem = malloc(sizeof(double)*k0->n); + + dcopy_(&mn, k0->a, &incr, mn_mem, &incr); + dcopy_(&m, mz->a, &incr, m_mem, &incr); + nnls_c(mn_mem, &m, &m, &n, m_mem, f0, &rnorm, + n_mem, m_mem2, index, &mode); + free(mn_mem); + free(m_mem); + free(m_mem2); + free(n_mem); +} + + +void CalcSim1D(MATRIX *k, double *f, double *baseline, MATRIX *Msim, + INDATA *d, HEAD *head){ + int incr = 1; + double alpha = 1.0, beta = 0.0; + double *deltaM; + double one_ = -1.0; + + Msim->a = malloc(sizeof(double)*k->m); + deltaM = malloc(sizeof(double)*k->m); + dgemv_("N", &k->m, &k->n, &alpha, k->a, &k->m, f, + &incr, &beta, Msim->a, &incr); + if (strcmp(head->TXtype, "T2inf")==0){ + baseline[0] = f[k->n-1]; + } + + dcopy_(&k->m, Msim->a, &incr, deltaM, &incr); + daxpy_(&k->m, &one_, d->m.a, &incr, deltaM, &incr); + head->sigma = dnrm2_(&k->m, deltaM, &incr)/sqrt(k->m); + free(deltaM); + +} + +void CalcSim(MATRIX *k1, MATRIX *k2, double *f, double *baseline, MATRIX *Msim, + INDATA *d, HEAD *head){ + int i, j, k=0, mn = k1->m*k2->n; + int mm = k1->n * k2->n; + MATRIX F, KF; + div_t temp; + double *deltaM; + double one_ = -1.0; + int incr = 1; + + F.m = k1->n; + F.n = k2->n; + F.a = malloc(sizeof(double)*F.m*F.n); + deltaM = malloc(sizeof(double)*mm); + for (i=0; iTXtype, "T1inf") == 0){ + j = 0; + for (i=0; ia, &incr, deltaM, &incr); + daxpy_(&mm, &one_, d->m.a, &incr, deltaM, &incr); + head->sigma = dnrm2_(&mm, deltaM, &incr)/sqrt(mm); + free(F.a); + free(deltaM); +} + +void Copyright(){ + int i; + printf("\n"); + for (i=0; i<70; i++) printf("-"); + printf("\n"); + printf("1D- and 2D-ILT program, MAY 2014\n"); + printf("Author: T. Ohkubo (Chiba University)\n\n"); + printf("The program can not be distributed without author's permission.\n"); + printf("Use of this program is limited to non-profit or evaluation purposes."); + + printf("\n"); + for (i=0; i<70; i++) printf("-"); + printf("\n"); + printf("\n"); + +} + +int main(int argc, char** argv){ + INDATA d; + MATRIX k0, k1, k2, s1z, s2z, u1z, u2z, v1z, v2z; + MATRIX k1z, k2z, mz_tmp1, mz_tmp2, mz, msim; + MATRIX x, y; /* T1とT2の離散化データ */ + double *f, *baseline; + HEAD head; + Copyright(); + + SetParameter(&head); + SetArgment(argc, argv, &head); + MakeKspace(head.TXmin, head.TXmax, head.TXnum, &x, head.TXspace); + MakeKspace(head.TYmin, head.TYmax, head.TYnum, &y, head.TYspace); + if (head.dimension == 2){ + SetDataSize(&head, &d); + ReadData(&head, &d); + Confirm(&head, &d); + printf("Making kernel...\n"); + MakeKernel(&d.u, &x, head.TXtype, &k1); + MakeKernel(&d.v, &y, head.TYtype, &k2); + printf("SVD compressing...\n"); + printf(" Kernel size (K1): %d*%d\n", k1.m, k1.n); + printf(" Kernel size (K2): %d*%d\n", k2.m, k2.n); + SvdCompress(&k1, &u1z, &s1z, &v1z, head.sfactor); + SvdCompress(&k2, &u2z, &s2z, &v2z, head.sfactor); + printf("Making compress kernel...\n"); + MultiplyMatrix(&s1z, &v1z, &k1z, "N", "N"); + MultiplyMatrix(&s2z, &v2z, &k2z, "N", "N"); + printf(" Kernel size (K1): %d*%d\n", k1z.m, k1z.n); + printf(" Kernel size (K2): %d*%d\n", k2z.m, k2z.n); + head.s1 = k1z.m; head.N1 = k1.m; + head.s2 = k2z.m; head.N2 = k2.m; + MultiplyMatrix(&u1z, &d.m, &mz_tmp1, "T", "N"); + MultiplyMatrix(&mz_tmp1, &u2z, &mz_tmp2, "N", "N"); + Lexicographic(&mz_tmp2, &mz); + Kronecker(&k1z, &k2z, &k0); + /* mz, k0以外のmallocしたメモリをrelease */ + free(s1z.a); free(s2z.a); + free(v1z.a); free(v2z.a); + free(u1z.a); free(u2z.a); + free(k1z.a); free(k2z.a); + free(mz_tmp1.a); free(mz_tmp2.a); + } + if (head.dimension == 1){ + ReadData1D(&head, &d); + Confirm(&head, &d); + printf("Making kernel...\n"); + MakeKernel(&d.u, &x, head.TXtype, &k1); + printf("Making compress kernel...\n"); + SvdCompress(&k1, &u1z, &s1z, &v1z, head.sfactor); + MultiplyMatrix(&s1z, &v1z, &k0, "N", "N"); + head.s1 = k0.m; head.N1 = k1.m; + k2.n = 1; + printf(" Kernel size (K1): %d*%d\n", k0.m, k0.n); + MultiplyMatrix(&u1z, &d.m, &mz, "T", "N"); + } + printf("Calculating initial F using NNLS... \n"); + f = malloc(sizeof(double)*k0.n); + baseline = malloc(sizeof(double)*k2.n); + + Calcf0(&k0, &mz, f); + + if (head.s_taget == 0) { + printf("Start BRD routine...\n"); + BRDAlphaLoop(&k0, &mz, &head, f, &k1, &k2, &d); + } + else { + printf("Start S-curve routine...\n"); + ScurveAlphaLoop(&k0, &mz, &head, f, &k1, &k2, &d); + } + printf("Output results...\n"); + if (head.dimension == 2){ + CalcSim(&k1, &k2, f, baseline, &msim, &d, &head); + OutputDistr(&head, &x, &y, f, baseline); + OutputDecay(&head, &d, &x, &y, &msim); + } + if (head.dimension == 1){ + CalcSim1D(&k1, f, baseline, &msim, &d, &head); + OutputDistr1D(&head, &x, f, baseline); + OutputDecay1D(&head, &d, &x, &msim); + } + Copyright(); + return 0; +} diff --git a/ilt.h b/ilt.h new file mode 100644 index 0000000..715e2cf --- /dev/null +++ b/ilt.h @@ -0,0 +1,90 @@ +#include +#include +#include +#include + +typedef struct{ + double newton_tol; + double alpha_tol; + double sfactor; + double alpha0; + double sigma0; + int confirm; + int alpha_loop_max; + int newton_loop_max; + int TXnum, TYnum; + int DXnum, DYnum; + int dimension; + double TXmin, TXmax; + double TYmin, TYmax; + int Atype; + char TXspace[10]; + char TYspace[10]; + char TXtype[10]; + char TYtype[10]; + char infile[1024]; + char outbase[1024]; + int newton_loop; + int alpha_loop; + /* s-curve用パラメータ */ + double s_taget; + double alpha_start; + double alpha_end; + int alpha_num; + double *alpha_list; + double *sigma_list; + + /* 計算結果用パラメータ */ + double alpha; + int s1, s2; + int N1, N2; + double sigma; + FILE* infp; +}HEAD; + +typedef struct{ + double *a; + int m, n; +}MATRIX; + +typedef struct{ + MATRIX u; /* tau1のベクトル */ + MATRIX v; /* tau2のベクトル */ + MATRIX m; /* 観測データの行列 u.m*v.m */ + double mmax; /* 観測データの最大値 */ +}INDATA; + +double ddot_(); +double dnrm2_(); +void dcopy_(); +void dgemm_(); +void dgesvd_(); +void dgesv_(); +void dgeev_(); +void dsysv_(); +void dgemv_(); +void nnls_c(); +void daxpy_(); +void dscal_(); +int idamax_(); +int SetArgment(); +void Transpose(); +void OutputDistr(); +void OutputDecay(); +void OutputDistr1D(); +void OutputDecay1D(); +void Error(); +void Confirm(); +void MakeKernel(); +void MultiplyMatrix(); +void Lexicographic(); +void ReadData(); +void ReadData1D(); +void SetDataSize(); +double BRDnewton_modify(); +void BRDcg(); +void BRDAlphaLoop(); +void ScurveAlphaLoop(); +void CalcSim(); +double CalcSigma(); +void OutputScurve(); diff --git a/input.c b/input.c new file mode 100644 index 0000000..1c7b26d --- /dev/null +++ b/input.c @@ -0,0 +1,140 @@ +#include "ilt.h" + +void SetDataSize(HEAD *head, INDATA *d){ + char buf[1024]; + char *tp; + + if((head->infp = fopen(head->infile, "r"))==NULL){ + printf("%s is not found.\n", head->infile); + exit(1); + } + + tp = fgets(buf, sizeof(buf), head->infp); + if (tp == NULL){ + printf("Read error. %s\n", head->infile); + exit(1); + } + rewind(head->infp); + + if (head->DXnum == 0 && head->DYnum == 0){ + /* スペース.を区切りに文字列を抽出 */ + tp = strtok(buf, "X:/Y:" ); + tp = strtok(NULL, "X:/Y:" ); + if (tp == NULL){ + printf("X Dimesion for input data error!\n"); + exit(1); + } + head->DXnum = atoi(tp); /* tau1軸がu */; + tp = strtok(NULL, "X:/Y:" ); + if (tp == NULL){ + printf("Y Dimesion for input data error!\n"); + exit(1); + } + head->DYnum = atoi(tp); /* tau2軸がv */ + } + if (head->DXnum == 0 || head->DYnum == 0){ + printf("XY Dimesion for input data error!\n"); + exit(1); + } + d->u.m = head->DXnum; + d->v.m = head->DYnum; + d->u.n = 1; + d->v.n = 1; + d->u.a = (double*)malloc(sizeof(double) * d->u.m ); + d->v.a = (double*)malloc(sizeof(double) * d->v.m); + d->m.a = (double*)malloc(sizeof(double) * d->u.m * d->v.m); + d->m.m = d->u.m; /* T1軸が行 */ + d->m.n = d->v.m; /* T2軸が列 */ +} + + +/* ####################################################################### */ +/* 1Dの測定データを読む */ +/* ####################################################################### */ +void ReadData1D(HEAD *head, INDATA *d){ + char buf[1024]; + int pnt = 0, col; + + if ((strcmp(head->infile, "-")==0)) head->infp = stdin; + else if((head->infp = fopen(head->infile, "r"))==NULL){ + printf("%s is not found.\n", head->infile); + exit(1); + } + if (head->DXnum == 0){ + while (fgets(buf, sizeof(buf), head->infp) != NULL){ + if(strstr(buf, "#") == NULL && buf[0] != '\n'){ + pnt++; + } + } + head->DXnum = pnt; + rewind(head->infp); + } + d->u.m = head->DXnum; + d->u.n = 1; + d->v.m = 0; + d->v.n = 0; + d->u.a = (double*)malloc(sizeof(double) * head->DXnum); + d->m.a = (double*)malloc(sizeof(double) * head->DXnum); + d->m.m = d->u.m; + d->m.n = 1; + head->s1 = head->DXnum; + head->s2 = 0; + + pnt = 0; + while (fgets(buf, sizeof(buf), head->infp) != NULL){ + if(strstr(buf, "#") == NULL && buf[0] != '\n' && buf[0] != '\r'){ + col = sscanf(buf, "%lf %lf", &d->u.a[pnt], &d->m.a[pnt]); + if (col != 2){ + printf("read error. line numeber=%d (%s)\n", pnt, head->infile); + exit(1); + } + /* d->u.a[pnt] = d->u.a[pnt] * 1e-3; /\* s->ms *\/ */ + + pnt++; + if (pnt == head->DXnum) break; + } + } + + fclose(head->infp); + + /* 最大値で規格化 */ + /* d->mmax = fabs(d->m.a[idamax_(&m, d->m.a, &incr)-1]); */ + /* alpha = 1.0/d->mmax; */ + /* dscal_(&m, &alpha, d->m.a, &incr); */ +} + +/* ####################################################################### */ +/* 測定データの読み込み */ +/* ####################################################################### */ +void ReadData(HEAD *head, INDATA *d){ + char buf[1024]; + int i=0, j=0, p, col; + double ta, tb, mi; + div_t temp; + + while (fgets(buf, sizeof(buf), head->infp) != NULL){ + if(strstr(buf, "#") == NULL && buf[0] != '\n' && buf[0] != '\r'){ + col = sscanf(buf, "%lf %lf %lf", &ta, &tb, &mi); + if (col != 3){ + printf("read error. line:%d (%s)\n", i, head->infile); + exit(1); + } + /* X軸一定でY軸変化のためM行列を転置 */ + temp = div(i, d->m.n); + p = temp.rem * d->m.m + temp.quot; + /* Mの並びがT1変化が先なら p=iにする */ + d->m.a[p] = mi; + if (i < d->v.m) d->v.a[i] = tb; + if ((i % d->v.m)==0) { + d->u.a[j] = ta; + j++; + } + i++; + } + } + + /* mn = d->m.m * d->m.n; */ + /* d->mmax = fabs(d->m.a[idamax_(&mn, d->m.a, &incr)-1]); */ + /* alpha = 1.0/d->mmax; */ + /* dscal_(&mn, &alpha, d->m.a, &incr); */ +} diff --git a/nnls.c b/nnls.c new file mode 100644 index 0000000..b9fb516 --- /dev/null +++ b/nnls.c @@ -0,0 +1,764 @@ +/* $Id: nnls.c,v 1.1 2002/02/21 22:56:45 tgkolda Exp $ */ +/* $Source: /space/CVS-APPSPACK//appspack-3/src/nnls.c,v $ */ + +/* Distributed with ASYNCHRONOUS PARALLEL PATTERN SEARCH (APPS) */ + +/* The routines in this file have been translated from Fortran to C by + f2c. Additional modifications have been made to remove the + dependencies on the f2c header file and library. The original + Fortran 77 code accompanies the SIAM Publications printing of + "Solving Least Squares Problems," by C. Lawson and R. Hanson and is + freely available at www.netlib.org/lawson-hanson/all. */ + +/* nnls.F -- translated by f2c (version 19970805). + You must link the resulting object file with the libraries: + -lf2c -lm (in that order) +*/ + +/* The next line was removed after the f2c translation */ +/* #include "f2c.h" */ + +/* The next lines were added after the f2c translation. Also swapped + abs for nnls_abs and max for nnls_max to avoid confusion with some + compilers. */ +#include +#include +#include +#define nnls_max(a,b) ((a) >= (b) ? (a) : (b)) +#define nnls_abs(x) ((x) >= 0 ? (x) : -(x)) +typedef int integer; +typedef double doublereal; + +/* The following subroutine was added after the f2c translation */ +double d_sign(double *a, double *b) +{ + double x; + x = (*a >= 0 ? *a : - *a); + return( *b >= 0 ? x : -x); +} + +/* Table of constant values */ + +static integer c__1 = 1; +static integer c__0 = 0; +static integer c__2 = 2; + +/* SUBROUTINE NNLS (A,MDA,M,N,B,X,RNORM,W,ZZ,INDEX,MODE) */ + +/* Algorithm NNLS: NONNEGATIVE LEAST SQUARES */ + +/* The original version of this code was developed by */ +/* Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory */ +/* 1973 JUN 15, and published in the book */ +/* "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. */ +/* Revised FEB 1995 to accompany reprinting of the book by SIAM. */ + +/* GIVEN AN M BY N MATRIX, A, AND AN M-VECTOR, B, COMPUTE AN */ +/* N-VECTOR, X, THAT SOLVES THE LEAST SQUARES PROBLEM */ + +/* A * X = B SUBJECT TO X .GE. 0 */ +/* ------------------------------------------------------------------ */ +/* Subroutine Arguments */ + +/* A(),MDA,M,N MDA IS THE FIRST DIMENSIONING PARAMETER FOR THE */ +/* ARRAY, A(). ON ENTRY A() CONTAINS THE M BY N */ +/* MATRIX, A. ON EXIT A() CONTAINS */ +/* THE PRODUCT MATRIX, Q*A , WHERE Q IS AN */ +/* M BY M ORTHOGONAL MATRIX GENERATED IMPLICITLY BY */ +/* THIS SUBROUTINE. */ +/* B() ON ENTRY B() CONTAINS THE M-VECTOR, B. ON EXIT B() CON- */ +/* TAINS Q*B. */ +/* X() ON ENTRY X() NEED NOT BE INITIALIZED. ON EXIT X() WILL */ +/* CONTAIN THE SOLUTION VECTOR. */ +/* RNORM ON EXIT RNORM CONTAINS THE EUCLIDEAN NORM OF THE */ +/* RESIDUAL VECTOR. */ +/* W() AN N-ARRAY OF WORKING SPACE. ON EXIT W() WILL CONTAIN */ +/* THE DUAL SOLUTION VECTOR. W WILL SATISFY W(I) = 0. */ +/* FOR ALL I IN SET P AND W(I) .LE. 0. FOR ALL I IN SET Z */ +/* ZZ() AN M-ARRAY OF WORKING SPACE. */ +/* INDEX() AN INTEGER WORKING ARRAY OF LENGTH AT LEAST N. */ +/* ON EXIT THE CONTENTS OF THIS ARRAY DEFINE THE SETS */ +/* P AND Z AS FOLLOWS.. */ + +/* INDEX(1) THRU INDEX(NSETP) = SET P. */ +/* INDEX(IZ1) THRU INDEX(IZ2) = SET Z. */ +/* IZ1 = NSETP + 1 = NPP1 */ +/* IZ2 = N */ +/* MODE THIS IS A SUCCESS-FAILURE FLAG WITH THE FOLLOWING */ +/* MEANINGS. */ +/* 1 THE SOLUTION HAS BEEN COMPUTED SUCCESSFULLY. */ +/* 2 THE DIMENSIONS OF THE PROBLEM ARE BAD. */ +/* EITHER M .LE. 0 OR N .LE. 0. */ +/* 3 ITERATION COUNT EXCEEDED. MORE THAN 3*N ITERATIONS. */ + +/* ------------------------------------------------------------------ */ +/* Subroutine */ int nnls_(a, mda, m, n, b, x, rnorm, w, zz, index, mode) + doublereal *a; + integer *mda, *m, *n; + doublereal *b, *x, *rnorm, *w, *zz; + integer *index, *mode; +{ + /* System generated locals */ + integer a_dim1, a_offset, i__1, i__2; + doublereal d__1, d__2; + + /* Builtin functions */ + /* The following lines were commented out after the f2c translation */ + /* double sqrt(); */ + /* integer s_wsfe(), do_fio(), e_wsfe(); */ + + /* Local variables */ + extern doublereal diff_(); + static integer iter; + static doublereal temp, wmax; + static integer i__, j, l; + static doublereal t, alpha, asave; + static integer itmax, izmax, nsetp; + extern /* Subroutine */ int g1_(); + static doublereal dummy, unorm, ztest, cc; + extern /* Subroutine */ int h12_(); + static integer ii, jj, ip; + static doublereal sm; + static integer iz, jz; + static doublereal up, ss; + static integer rtnkey, iz1, iz2, npp1; + + /* Fortran I/O blocks */ + /* The following line was commented out after the f2c translation */ + /* static cilist io___22 = { 0, 6, 0, "(/a)", 0 }; */ + + + /* ------------------------------------------------------------------ + */ + /* integer INDEX(N) */ + /* double precision A(MDA,N), B(M), W(N), X(N), ZZ(M) */ + /* ------------------------------------------------------------------ + */ + /* Parameter adjustments */ + a_dim1 = *mda; + a_offset = a_dim1 + 1; + a -= a_offset; + --b; + --x; + --w; + --zz; + --index; + + /* Function Body */ + *mode = 1; + if (*m <= 0 || *n <= 0) { + *mode = 2; + return 0; + } + iter = 0; + itmax = *n * 3; + + /* INITIALIZE THE ARRAYS INDEX() AND X(). */ + + i__1 = *n; + for (i__ = 1; i__ <= i__1; ++i__) { + x[i__] = 0.; + /* L20: */ + index[i__] = i__; + } + + iz2 = *n; + iz1 = 1; + nsetp = 0; + npp1 = 1; + /* ****** MAIN LOOP BEGINS HERE ****** */ + L30: + /* QUIT IF ALL COEFFICIENTS ARE ALREADY IN THE SOLUTION. + */ + /* OR IF M COLS OF A HAVE BEEN TRIANGULARIZED. */ + + if (iz1 > iz2 || nsetp >= *m) { + goto L350; + } + + /* COMPUTE COMPONENTS OF THE DUAL (NEGATIVE GRADIENT) VECTOR W(). + */ + + i__1 = iz2; + for (iz = iz1; iz <= i__1; ++iz) { + j = index[iz]; + sm = 0.; + i__2 = *m; + for (l = npp1; l <= i__2; ++l) { + /* L40: */ + sm += a[l + j * a_dim1] * b[l]; + } + w[j] = sm; + /* L50: */ + } + /* FIND LARGEST POSITIVE W(J). */ + L60: + wmax = 0.; + i__1 = iz2; + for (iz = iz1; iz <= i__1; ++iz) { + j = index[iz]; + if (w[j] > wmax) { + wmax = w[j]; + izmax = iz; + } + /* L70: */ + } + + /* IF WMAX .LE. 0. GO TO TERMINATION. */ + /* THIS INDICATES SATISFACTION OF THE KUHN-TUCKER CONDITIONS. + */ + + if (wmax <= 0.) { + goto L350; + } + iz = izmax; + j = index[iz]; + + /* THE SIGN OF W(J) IS OK FOR J TO BE MOVED TO SET P. */ + /* BEGIN THE TRANSFORMATION AND CHECK NEW DIAGONAL ELEMENT TO AVOID */ + /* NEAR LINEAR DEPENDENCE. */ + + asave = a[npp1 + j * a_dim1]; + i__1 = npp1 + 1; + h12_(&c__1, &npp1, &i__1, m, &a[j * a_dim1 + 1], &c__1, &up, &dummy, & + c__1, &c__1, &c__0); + unorm = 0.; + if (nsetp != 0) { + i__1 = nsetp; + for (l = 1; l <= i__1; ++l) { + /* L90: */ + /* Computing 2nd power */ + d__1 = a[l + j * a_dim1]; + unorm += d__1 * d__1; + } + } + unorm = sqrt(unorm); + d__2 = unorm + (d__1 = a[npp1 + j * a_dim1], nnls_abs(d__1)) * .01; + if (diff_(&d__2, &unorm) > 0.) { + + /* COL J IS SUFFICIENTLY INDEPENDENT. COPY B INTO ZZ, UPDATE Z + Z */ + /* AND SOLVE FOR ZTEST ( = PROPOSED NEW VALUE FOR X(J) ). */ + + i__1 = *m; + for (l = 1; l <= i__1; ++l) { + /* L120: */ + zz[l] = b[l]; + } + i__1 = npp1 + 1; + h12_(&c__2, &npp1, &i__1, m, &a[j * a_dim1 + 1], &c__1, &up, &zz[1], & + c__1, &c__1, &c__1); + ztest = zz[npp1] / a[npp1 + j * a_dim1]; + + /* SEE IF ZTEST IS POSITIVE */ + + if (ztest > 0.) { + goto L140; + } + } + + /* REJECT J AS A CANDIDATE TO BE MOVED FROM SET Z TO SET P. */ + /* RESTORE A(NPP1,J), SET W(J)=0., AND LOOP BACK TO TEST DUAL */ + /* COEFFS AGAIN. */ + + a[npp1 + j * a_dim1] = asave; + w[j] = 0.; + goto L60; + + /* THE INDEX J=INDEX(IZ) HAS BEEN SELECTED TO BE MOVED FROM */ + /* SET Z TO SET P. UPDATE B, UPDATE INDICES, APPLY HOUSEHOLDER */ + /* TRANSFORMATIONS TO COLS IN NEW SET Z, ZERO SUBDIAGONAL ELTS IN */ + /* COL J, SET W(J)=0. */ + + L140: + i__1 = *m; + for (l = 1; l <= i__1; ++l) { + /* L150: */ + b[l] = zz[l]; + } + + index[iz] = index[iz1]; + index[iz1] = j; + ++iz1; + nsetp = npp1; + ++npp1; + + if (iz1 <= iz2) { + i__1 = iz2; + for (jz = iz1; jz <= i__1; ++jz) { + jj = index[jz]; + h12_(&c__2, &nsetp, &npp1, m, &a[j * a_dim1 + 1], &c__1, &up, &a[ + jj * a_dim1 + 1], &c__1, mda, &c__1); + /* L160: */ + } + } + + if (nsetp != *m) { + i__1 = *m; + for (l = npp1; l <= i__1; ++l) { + /* L180: */ + a[l + j * a_dim1] = 0.; + } + } + + w[j] = 0.; + /* SOLVE THE TRIANGULAR SYSTEM. */ + /* STORE THE SOLUTION TEMPORARILY IN ZZ(). + */ + rtnkey = 1; + goto L400; + L200: + + /* ****** SECONDARY LOOP BEGINS HERE ****** */ + + /* ITERATION COUNTER. */ + + L210: + ++iter; + if (iter > itmax) { + *mode = 3; + /* The following lines were replaced after the f2c translation */ + /* s_wsfe(&io___22); */ + /* do_fio(&c__1, " NNLS quitting on iteration count.", 34L); */ + /* e_wsfe(); */ + fprintf(stdout, "\n NNLS quitting on iteration count.\n"); + fflush(stdout); + goto L350; + } + + /* SEE IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE. */ + /* IF NOT COMPUTE ALPHA. */ + + alpha = 2.; + i__1 = nsetp; + for (ip = 1; ip <= i__1; ++ip) { + l = index[ip]; + if (zz[ip] <= 0.) { + t = -x[l] / (zz[ip] - x[l]); + if (alpha > t) { + alpha = t; + jj = ip; + } + } + /* L240: */ + } + + /* IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE THEN ALPHA WILL */ + /* STILL = 2. IF SO EXIT FROM SECONDARY LOOP TO MAIN LOOP. */ + + if (alpha == 2.) { + goto L330; + } + + /* OTHERWISE USE ALPHA WHICH WILL BE BETWEEN 0. AND 1. TO */ + /* INTERPOLATE BETWEEN THE OLD X AND THE NEW ZZ. */ + + i__1 = nsetp; + for (ip = 1; ip <= i__1; ++ip) { + l = index[ip]; + x[l] += alpha * (zz[ip] - x[l]); + /* L250: */ + } + + /* MODIFY A AND B AND THE INDEX ARRAYS TO MOVE COEFFICIENT I */ + /* FROM SET P TO SET Z. */ + + i__ = index[jj]; + L260: + x[i__] = 0.; + + if (jj != nsetp) { + ++jj; + i__1 = nsetp; + for (j = jj; j <= i__1; ++j) { + ii = index[j]; + index[j - 1] = ii; + g1_(&a[j - 1 + ii * a_dim1], &a[j + ii * a_dim1], &cc, &ss, &a[j + - 1 + ii * a_dim1]); + a[j + ii * a_dim1] = 0.; + i__2 = *n; + for (l = 1; l <= i__2; ++l) { + if (l != ii) { + + /* Apply procedure G2 (CC,SS,A(J-1,L),A(J, + L)) */ + + temp = a[j - 1 + l * a_dim1]; + a[j - 1 + l * a_dim1] = cc * temp + ss * a[j + l * a_dim1] + ; + a[j + l * a_dim1] = -ss * temp + cc * a[j + l * a_dim1]; + } + /* L270: */ + } + + /* Apply procedure G2 (CC,SS,B(J-1),B(J)) */ + + temp = b[j - 1]; + b[j - 1] = cc * temp + ss * b[j]; + b[j] = -ss * temp + cc * b[j]; + /* L280: */ + } + } + + npp1 = nsetp; + --nsetp; + --iz1; + index[iz1] = i__; + + /* SEE IF THE REMAINING COEFFS IN SET P ARE FEASIBLE. THEY SHOULD + */ + /* BE BECAUSE OF THE WAY ALPHA WAS DETERMINED. */ + /* IF ANY ARE INFEASIBLE IT IS DUE TO ROUND-OFF ERROR. ANY */ + /* THAT ARE NONPOSITIVE WILL BE SET TO ZERO */ + /* AND MOVED FROM SET P TO SET Z. */ + + i__1 = nsetp; + for (jj = 1; jj <= i__1; ++jj) { + i__ = index[jj]; + if (x[i__] <= 0.) { + goto L260; + } + /* L300: */ + } + + /* COPY B( ) INTO ZZ( ). THEN SOLVE AGAIN AND LOOP BACK. */ + + i__1 = *m; + for (i__ = 1; i__ <= i__1; ++i__) { + /* L310: */ + zz[i__] = b[i__]; + } + rtnkey = 2; + goto L400; + L320: + goto L210; + /* ****** END OF SECONDARY LOOP ****** */ + + L330: + i__1 = nsetp; + for (ip = 1; ip <= i__1; ++ip) { + i__ = index[ip]; + /* L340: */ + x[i__] = zz[ip]; + } + /* ALL NEW COEFFS ARE POSITIVE. LOOP BACK TO BEGINNING. */ + goto L30; + + /* ****** END OF MAIN LOOP ****** */ + + /* COME TO HERE FOR TERMINATION. */ + /* COMPUTE THE NORM OF THE FINAL RESIDUAL VECTOR. */ + + L350: + sm = 0.; + if (npp1 <= *m) { + i__1 = *m; + for (i__ = npp1; i__ <= i__1; ++i__) { + /* L360: */ + /* Computing 2nd power */ + d__1 = b[i__]; + sm += d__1 * d__1; + } + } else { + i__1 = *n; + for (j = 1; j <= i__1; ++j) { + /* L380: */ + w[j] = 0.; + } + } + *rnorm = sqrt(sm); + return 0; + + /* THE FOLLOWING BLOCK OF CODE IS USED AS AN INTERNAL SUBROUTINE */ + /* TO SOLVE THE TRIANGULAR SYSTEM, PUTTING THE SOLUTION IN ZZ(). */ + + L400: + i__1 = nsetp; + for (l = 1; l <= i__1; ++l) { + ip = nsetp + 1 - l; + if (l != 1) { + i__2 = ip; + for (ii = 1; ii <= i__2; ++ii) { + zz[ii] -= a[ii + jj * a_dim1] * zz[ip + 1]; + /* L410: */ + } + } + jj = index[ip]; + zz[ip] /= a[ip + jj * a_dim1]; + /* L430: */ + } + switch ((int)rtnkey) { + case 1: goto L200; + case 2: goto L320; + } + + /* The next line was added after the f2c translation to keep + compilers from complaining about a void return from a non-void + function. */ + return 0; + +} /* nnls_ */ + +/* Subroutine */ int g1_(a, b, cterm, sterm, sig) + doublereal *a, *b, *cterm, *sterm, *sig; +{ + /* System generated locals */ + doublereal d__1; + + /* Builtin functions */ + /* The following line was commented out after the f2c translation */ + /* double sqrt(), d_sign(); */ + + /* Local variables */ + static doublereal xr, yr; + + + /* COMPUTE ORTHOGONAL ROTATION MATRIX.. */ + + /* The original version of this code was developed by */ + /* Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory + */ + /* 1973 JUN 12, and published in the book */ + /* "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. */ + /* Revised FEB 1995 to accompany reprinting of the book by SIAM. */ + + /* COMPUTE.. MATRIX (C, S) SO THAT (C, S)(A) = (SQRT(A**2+B**2)) */ + /* (-S,C) (-S,C)(B) ( 0 ) */ + /* COMPUTE SIG = SQRT(A**2+B**2) */ + /* SIG IS COMPUTED LAST TO ALLOW FOR THE POSSIBILITY THAT */ + /* SIG MAY BE IN THE SAME LOCATION AS A OR B . */ + /* ------------------------------------------------------------------ + */ + /* ------------------------------------------------------------------ + */ + if (nnls_abs(*a) > nnls_abs(*b)) { + xr = *b / *a; + /* Computing 2nd power */ + d__1 = xr; + yr = sqrt(d__1 * d__1 + 1.); + d__1 = 1. / yr; + *cterm = d_sign(&d__1, a); + *sterm = *cterm * xr; + *sig = nnls_abs(*a) * yr; + return 0; + } + if (*b != 0.) { + xr = *a / *b; + /* Computing 2nd power */ + d__1 = xr; + yr = sqrt(d__1 * d__1 + 1.); + d__1 = 1. / yr; + *sterm = d_sign(&d__1, b); + *cterm = *sterm * xr; + *sig = nnls_abs(*b) * yr; + return 0; + } + *sig = 0.; + *cterm = 0.; + *sterm = 1.; + return 0; +} /* g1_ */ + +/* SUBROUTINE H12 (MODE,LPIVOT,L1,M,U,IUE,UP,C,ICE,ICV,NCV) */ + +/* CONSTRUCTION AND/OR APPLICATION OF A SINGLE */ +/* HOUSEHOLDER TRANSFORMATION.. Q = I + U*(U**T)/B */ + +/* The original version of this code was developed by */ +/* Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory */ +/* 1973 JUN 12, and published in the book */ +/* "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. */ +/* Revised FEB 1995 to accompany reprinting of the book by SIAM. */ +/* ------------------------------------------------------------------ */ +/* Subroutine Arguments */ + +/* MODE = 1 OR 2 Selects Algorithm H1 to construct and apply a */ +/* Householder transformation, or Algorithm H2 to apply a */ +/* previously constructed transformation. */ +/* LPIVOT IS THE INDEX OF THE PIVOT ELEMENT. */ +/* L1,M IF L1 .LE. M THE TRANSFORMATION WILL BE CONSTRUCTED TO */ +/* ZERO ELEMENTS INDEXED FROM L1 THROUGH M. IF L1 GT. M */ +/* THE SUBROUTINE DOES AN IDENTITY TRANSFORMATION. */ +/* U(),IUE,UP On entry with MODE = 1, U() contains the pivot */ +/* vector. IUE is the storage increment between elements. */ +/* On exit when MODE = 1, U() and UP contain quantities */ +/* defining the vector U of the Householder transformation. */ +/* on entry with MODE = 2, U() and UP should contain */ +/* quantities previously computed with MODE = 1. These will */ +/* not be modified during the entry with MODE = 2. */ +/* C() ON ENTRY with MODE = 1 or 2, C() CONTAINS A MATRIX WHICH */ +/* WILL BE REGARDED AS A SET OF VECTORS TO WHICH THE */ +/* HOUSEHOLDER TRANSFORMATION IS TO BE APPLIED. */ +/* ON EXIT C() CONTAINS THE SET OF TRANSFORMED VECTORS. */ +/* ICE STORAGE INCREMENT BETWEEN ELEMENTS OF VECTORS IN C(). */ +/* ICV STORAGE INCREMENT BETWEEN VECTORS IN C(). */ +/* NCV NUMBER OF VECTORS IN C() TO BE TRANSFORMED. IF NCV .LE. 0 */ +/* NO OPERATIONS WILL BE DONE ON C(). */ +/* ------------------------------------------------------------------ */ +/* Subroutine */ int h12_(mode, lpivot, l1, m, u, iue, up, c__, ice, icv, ncv) + integer *mode, *lpivot, *l1, *m; + doublereal *u; + integer *iue; + doublereal *up, *c__; + integer *ice, *icv, *ncv; +{ + /* System generated locals */ + integer u_dim1, u_offset, i__1, i__2; + doublereal d__1, d__2; + + /* Builtin functions */ + /* The following line was commented out after the f2c translation */ + /* double sqrt(); */ + + /* Local variables */ + static integer incr; + static doublereal b; + static integer i__, j; + static doublereal clinv; + static integer i2, i3, i4; + static doublereal cl, sm; + + /* ------------------------------------------------------------------ + */ + /* double precision U(IUE,M) */ + /* ------------------------------------------------------------------ + */ + /* Parameter adjustments */ + u_dim1 = *iue; + u_offset = u_dim1 + 1; + u -= u_offset; + --c__; + + /* Function Body */ + if (0 >= *lpivot || *lpivot >= *l1 || *l1 > *m) { + return 0; + } + cl = (d__1 = u[*lpivot * u_dim1 + 1], nnls_abs(d__1)); + if (*mode == 2) { + goto L60; + } + /* ****** CONSTRUCT THE TRANSFORMATION. ****** + */ + i__1 = *m; + for (j = *l1; j <= i__1; ++j) { + /* L10: */ + /* Computing MAX */ + d__2 = (d__1 = u[j * u_dim1 + 1], nnls_abs(d__1)); + cl = nnls_max(d__2,cl); + } + if (cl <= 0.) { + goto L130; + } else { + goto L20; + } + L20: + clinv = 1. / cl; + /* Computing 2nd power */ + d__1 = u[*lpivot * u_dim1 + 1] * clinv; + sm = d__1 * d__1; + i__1 = *m; + for (j = *l1; j <= i__1; ++j) { + /* L30: */ + /* Computing 2nd power */ + d__1 = u[j * u_dim1 + 1] * clinv; + sm += d__1 * d__1; + } + cl *= sqrt(sm); + if (u[*lpivot * u_dim1 + 1] <= 0.) { + goto L50; + } else { + goto L40; + } + L40: + cl = -cl; + L50: + *up = u[*lpivot * u_dim1 + 1] - cl; + u[*lpivot * u_dim1 + 1] = cl; + goto L70; + /* ****** APPLY THE TRANSFORMATION I+U*(U**T)/B TO C. ****** + */ + + L60: + if (cl <= 0.) { + goto L130; + } else { + goto L70; + } + L70: + if (*ncv <= 0) { + return 0; + } + b = *up * u[*lpivot * u_dim1 + 1]; + /* B MUST BE NONPOSITIVE HERE. IF B = 0., RETURN. + */ + + if (b >= 0.) { + goto L130; + } else { + goto L80; + } + L80: + b = 1. / b; + i2 = 1 - *icv + *ice * (*lpivot - 1); + incr = *ice * (*l1 - *lpivot); + i__1 = *ncv; + for (j = 1; j <= i__1; ++j) { + i2 += *icv; + i3 = i2 + incr; + i4 = i3; + sm = c__[i2] * *up; + i__2 = *m; + for (i__ = *l1; i__ <= i__2; ++i__) { + sm += c__[i3] * u[i__ * u_dim1 + 1]; + /* L90: */ + i3 += *ice; + } + if (sm != 0.) { + goto L100; + } else { + goto L120; + } + L100: + sm *= b; + c__[i2] += sm * *up; + i__2 = *m; + for (i__ = *l1; i__ <= i__2; ++i__) { + c__[i4] += sm * u[i__ * u_dim1 + 1]; + /* L110: */ + i4 += *ice; + } + L120: + ; + } + L130: + return 0; +} /* h12_ */ + +doublereal diff_(x, y) + doublereal *x, *y; +{ + /* System generated locals */ + doublereal ret_val; + + + /* Function used in tests that depend on machine precision. */ + + /* The original version of this code was developed by */ + /* Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory + */ + /* 1973 JUN 7, and published in the book */ + /* "SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974. */ + /* Revised FEB 1995 to accompany reprinting of the book by SIAM. */ + + ret_val = *x - *y; + return ret_val; +} /* diff_ */ + + +/* The following subroutine was added after the f2c translation */ +int nnls_c(double* a, const int* mda, const int* m, const int* n, double* b, + double* x, double* rnorm, double* w, double* zz, int* index, + int* mode) +{ + return (nnls_(a, mda, m, n, b, x, rnorm, w, zz, index, mode)); +} diff --git a/output.c b/output.c new file mode 100644 index 0000000..ce35dc5 --- /dev/null +++ b/output.c @@ -0,0 +1,208 @@ +#include "ilt.h" + +/* ####################################################################### */ +/* Make header information */ +/* ####################################################################### */ +void OutputHeader(FILE *fp, HEAD *head){ + fprintf(fp, "# -------------- solving condition -------------- #\n"); + fprintf(fp, "# sfactor: %e\n", head->sfactor); + fprintf(fp, "# initial alpha: %g\n", head->alpha0); + fprintf(fp, "# alpha max loop: %d\n", head->alpha_loop_max); + fprintf(fp, "# alpha tol: %g\n", head->alpha_tol); + fprintf(fp, "# Newton max loop: %d\n", head->newton_loop_max); + fprintf(fp, "# Newton tol: %e\n", head->newton_tol); + fprintf(fp, "# A_opt method: %d\n", head->Atype); + fprintf(fp, "# --------------- TX information ---------------- #\n"); + fprintf(fp, "# Data space: %s\n", head->TXspace); + fprintf(fp, "# Data type: %s\n", head->TXtype); + fprintf(fp, "# (TXmin, TXmax, TXpoints)="); + fprintf(fp, "(%g, %g, %d)\n", head->TXmin, head->TXmax, head->TXnum); + if (head->dimension == 2){ + fprintf(fp, "# ---------------- TY information --------------- #\n"); + fprintf(fp, "# Data space: %s\n", head->TYspace); + fprintf(fp, "# Data type: %s\n", head->TYtype); + fprintf(fp, "# (TYmin, TYmax, TYpoints)="); + fprintf(fp, "(%g, %g, %d)\n", head->TYmin, head->TYmax, head->TYnum); + } + fprintf(fp, "# ------------------ input data ----------------- #\n"); + fprintf(fp, "# Input file: %s\n", head->infile); + fprintf(fp, "# Output basename: %s\n", head->outbase); + fprintf(fp, "# ---------------- Compress data ---------------- #\n"); + fprintf(fp, "# K1 original : %d\n", head->N1); + fprintf(fp, "# K1 compress : %d\n", head->s1); + if (head->dimension == 2){ + fprintf(fp, "# K2 original : %d\n", head->N2); + fprintf(fp, "# K2 compress : %d\n", head->s2); + } + fprintf(fp, "# ------------- Calculation results ------------- #\n"); + fprintf(fp, "# Alpha : %e\n", head->alpha); + fprintf(fp, "# Sigma : %e\n", head->sigma); + fprintf(fp, "# ----------------------------------------- #\n"); + fprintf(fp, "\n"); +} + +/* ####################################################################### */ +/* output for distribution data*/ +/* ####################################################################### */ +void OutputDistr(HEAD *head, MATRIX *x, MATRIX *y, double *f, double *baseline){ + int i, j, k=0; + char outfile[2048]; + FILE *fp; + double dmax; + int incr = 1; + int mn = x->m * y->m; + + printf(" Distribution data: %s.distr\n", head->outbase); + if ((strcmp(head->outbase, "-")==0)) fp = stdout; + else { + sprintf(outfile, "%s.distr", head->outbase); + fp = fopen(outfile, "w"); + } + OutputHeader(fp, head); + dmax = fabs(f[idamax_(&mn, f, &incr)-1]); + fprintf(fp, "# Distribution normalized coeff.: %g\n", dmax); + fprintf(fp, "# %13s %15s %15s", "TX", "TY", "I\n"); + for (i=0; im; i++){ + for (j=0; jm; j++){ + fprintf(fp, "%15.8e %15.8e %15.8e\n", x->a[i], y->a[j], f[k]/dmax); + k++; + } + fprintf(fp, "\n"); + } + fprintf(fp, "\n\n\n"); + fprintf(fp, "# Baseline list\n"); + for (i=0; iTYnum; i++){ + fprintf(fp, "# %15.8e\n", baseline[i]/dmax); + } + fclose(fp); +} + +/* ####################################################################### */ +/* output for decay data */ +/* ####################################################################### */ +void OutputDecay(HEAD *head, INDATA *d, + MATRIX *x, MATRIX *y, MATRIX *m){ + int i, j, p=0; + double t1, t2; + char outfile[2048]; + MATRIX Iexp, Isim; + FILE *fp; + + printf(" Decay data: %s%s\n", head->outbase, ".decay"); + Transpose(&d->m, &Iexp); + Transpose(m, &Isim); + if ((strcmp(head->outbase, "-")==0)) fp = stdout; + else { + sprintf(outfile, "%s.decay", head->outbase); + fp = fopen(outfile, "w"); + } + OutputHeader(fp, head); + fprintf(fp, "# Signal normalized coeff.: %g\n", d->mmax); + fprintf(fp, "# %13s %15s %15s %15s %15s\n", + "t1", "t2", "Iexp", "Isim", "Iexp - Isim"); + for (i=0; iu.m; i++){ + for (j=0; jv.m; j++){ + t1 = d->u.a[i]; + t2 = d->v.a[j]; + fprintf(fp, "%15.8e %15.8e %15.8e %15.8e %15.8e\n", + t1, t2, Iexp.a[p], Isim.a[p], Iexp.a[p] - Isim.a[p]); + p++; + } + fprintf(fp, "\n\n"); + } + fclose(fp); + if (head->s_taget != 0) OutputScurve(head); +} + +/* ####################################################################### */ +/* output for distribution data */ +/* ####################################################################### */ +void OutputDistr1D(HEAD *head, MATRIX *x, double *f, double *baseline){ + int i; + char outfile[2048]; + FILE *fp; + double dmax; + int incr = 1; + int mn = x->m; + + printf(" Distribution data: %s.distr\n", head->outbase); + if ((strcmp(head->outbase, "-")==0)) fp = stdout; + else { + sprintf(outfile, "%s%s", head->outbase, ".distr"); + fp = fopen(outfile, "w"); + } + OutputHeader(fp, head); + dmax = fabs(f[idamax_(&mn, f, &incr)-1]); + fprintf(fp, "# Distribution normalized coeff.: %g\n", dmax); + fprintf(fp, "# %13s %15s \n", "TX", "I"); + for (i=0; im; i++){ + fprintf(fp, "%15.8e %15.8e\n", x->a[i], f[i]/dmax); + } + fprintf(fp, "\n\n\n"); + fprintf(fp, "# Baseline list\n"); + fprintf(fp, "# %15.8e\n", baseline[0]/dmax); + if ((strcmp(head->outbase, "-")!=0)) fclose(fp); +} + +/* ####################################################################### */ +/* output for decay data */ +/* ####################################################################### */ +void OutputDecay1D(HEAD *head, INDATA *d, MATRIX *x, MATRIX *m){ + int i; + double t1; + char outfile[2048]; + FILE *fp; + + printf(" Decay data: %s.decay\n", head->outbase); + if ((strcmp(head->outbase, "-")==0)) fp = stdout; + else { + sprintf(outfile, "%s.decay", head->outbase); + fp = fopen(outfile, "w"); + } + OutputHeader(fp, head); + fprintf(fp, "# Signal normalized coeff.: %g\n", d->mmax); + fprintf(fp, "# %13s %15s %15s %15s\n", + "t1", "Iexp", "Isim", "Iexp - Isim"); + for (i=0; iu.m; i++){ + t1 = d->u.a[i]; + fprintf(fp, "%15.8e %15.8e %15.8e %15.8e\n", + t1, d->m.a[i], m->a[i], d->m.a[i] - m->a[i]); + } + fclose(fp); + fflush(fp); +} + +/* ####################################################################### */ +/* output for Scurve */ +/* ####################################################################### */ +void OutputScurve(HEAD *head){ + int i, sw=0; + double delta; + char outfile[2048]; + FILE *fp; + + printf(" Scurve data: %s.scurve\n", head->outbase); + if ((strcmp(head->outbase, "-")==0)) fp = stdout; + else { + sprintf(outfile, "%s.scurve", head->outbase); + fp = fopen(outfile, "w"); + } + fprintf(fp, "# %1s %15s %15s %15s %15s %15s\n", + "n", "alpha", "sigma", "log10(alpha)", "log10(sigma)", "grad"); + for (i=0; ialpha_num; i++){ + if (i!=0 && (head->alpha_list[i] < head->alpha_list[i-1]) && sw==0){ + fprintf(fp, "\n\n"); + sw = 1; + } + fprintf(fp, "%3d %15.8e %15.8e %15.8e %15.8e", + i+1, head->alpha_list[i], head->sigma_list[i], + log10(head->alpha_list[i]), log10(head->sigma_list[i])); + if (i!=(head->alpha_num-1)) { + delta = (log10(head->sigma_list[i+1]) - log10(head->sigma_list[i]))/ + (log10(head->alpha_list[i+1]) - log10(head->alpha_list[i])); + fprintf(fp, " %15.8e\n", delta); + } + else fprintf(fp, "\n"); + } +} +