/* armaest22.f -- translated by f2c (version 19991025). You must link the resulting object file with the libraries: -lf2c -lm (in that order) */ #include "f2c.h" /* Common Block Declarations */ struct { integer nar, nma, arlst[100] /* was [4][25] */, malst[100] /* was [4][25] */; } dat2_; #define dat2_1 dat2_ struct { char mdltst[5]; } diag_; #define diag_1 diag_ /* Table of constant values */ static integer c__9 = 9; static integer c__1 = 1; static integer c__3 = 3; static integer c__2 = 2; static integer c__5 = 5; /* Main program */ MAIN__() { /* Format strings */ static char fmt_6[] = "(a80)"; static char fmt_77[] = "(i4)"; static char fmt_18[] = "(i16)"; static char fmt_691[] = "(i6,i6)"; static char fmt_54[] = "(/)"; static char fmt_431[] = "(2i4,4x,e21.14)"; static char fmt_66[] = "(i3)"; static char fmt_2[] = "(2i4)"; static char fmt_21[] = "(2i4/2i4)"; static char fmt_43[] = "(2i4,4x,e14.7,4x,e14.7)"; static char fmt_421[] = "(\002AR parameter estimates\002)"; static char fmt_42[] = "(\002MA parameter estimates\002)"; static char fmt_63[] = "(2i4,4x,e13.6,2x,e13.6,4x,e13.6,2x,e13.6)"; /* System generated locals */ address a__1[2]; integer i__1[2], i__2, i__3, i__4; doublereal d__1; olist o__1; cllist cl__1; alist al__1; /* Builtin functions */ integer s_rsle(), do_lio(), e_rsle(), s_wsle(), e_wsle(), s_cmp(); /* Subroutine */ int s_cat(); integer f_open(), s_rsfe(), do_fio(), e_rsfe(), f_clos(), f_rew(), s_wsfe( ), e_wsfe(); /* Subroutine */ int s_copy(); double sqrt(); /* Local variables */ static doublereal sdma[100] /* was [25][4] */, sdar[100] /* was [25][4] */; static integer ncol; static char path[80]; static integer nset, nspl, ncol1; static doublereal gtol1; static integer pset1, pset2, pset3, c__[4] /* was [2][2] */, nrow1; static doublereal h__[262144] /* was [512][512] */; static integer i__, j, k, l; static char label[30]; static doublereal p[25]; static integer q[4] /* was [2][2] */; static char maflg[20]; static doublereal x[262144] /* was [512][512] */; static integer iseed; static char adnse[5], sdtbl[20], input[20], label1[30]; static integer i8; static char label3[30]; static doublereal xyvar; extern /* Subroutine */ int dfinv1_(); static integer malst0[100] /* was [4][25] */; static doublereal ma[5000] /* was [25][4][50] */, ar[5000] /* was [25][4] [50] */; static integer sd[50]; static shortint xr[262144] /* was [512][512] */; static char ftrgam[30]; static doublereal xy[304704] /* was [552][552] */; static char ftrflg[10]; static doublereal ma0[1600] /* was [25][64] */, ma1[25], ar1[1600] /* was [25][64] */, ar2[25]; static char string[80]; static integer il2, ik2, ik3, il3; extern /* Subroutine */ int imgmdl2_(); static doublereal xy1[304704] /* was [552][552] */; static integer ndb; extern /* Subroutine */ int imprsp2_(); static integer nmc; static doublereal uma[100] /* was [25][4] */; static integer cnr[32] /* was [2][16] */; static doublereal nrf, uar[100] /* was [25][4] */; static char nse[5], labelim[30]; static integer nsd; static char est[5], blurflg[20]; static doublereal uxy; static integer nma0; /* Fortran I/O blocks */ static cilist io___1 = { 0, 5, 0, 0, 0 }; static cilist io___2 = { 0, 6, 0, 0, 0 }; static cilist io___3 = { 0, 5, 0, 0, 0 }; static cilist io___6 = { 0, 6, 0, 0, 0 }; static cilist io___7 = { 0, 5, 0, 0, 0 }; static cilist io___11 = { 0, 6, 0, 0, 0 }; static cilist io___12 = { 0, 5, 0, 0, 0 }; static cilist io___14 = { 0, 6, 0, 0, 0 }; static cilist io___15 = { 0, 5, 0, 0, 0 }; static cilist io___16 = { 0, 6, 0, 0, 0 }; static cilist io___17 = { 0, 5, 0, 0, 0 }; static cilist io___19 = { 0, 6, 0, 0, 0 }; static cilist io___20 = { 0, 6, 0, 0, 0 }; static cilist io___21 = { 0, 5, 0, 0, 0 }; static cilist io___23 = { 0, 6, 0, 0, 0 }; static cilist io___24 = { 0, 6, 0, 0, 0 }; static cilist io___25 = { 0, 5, 0, 0, 0 }; static cilist io___27 = { 0, 6, 0, 0, 0 }; static cilist io___29 = { 0, 3, 0, fmt_6, 0 }; static cilist io___31 = { 0, 3, 0, fmt_77, 0 }; static cilist io___33 = { 0, 3, 0, fmt_18, 0 }; static cilist io___36 = { 0, 5, 0, 0, 0 }; static cilist io___39 = { 0, 6, 0, 0, 0 }; static cilist io___40 = { 0, 3, 0, fmt_6, 0 }; static cilist io___41 = { 0, 3, 0, fmt_691, 0 }; static cilist io___43 = { 0, 3, 0, fmt_54, 0 }; static cilist io___45 = { 0, 3, 0, fmt_431, 0 }; static cilist io___48 = { 0, 6, 0, fmt_431, 0 }; static cilist io___49 = { 0, 3, 0, fmt_54, 0 }; static cilist io___50 = { 0, 6, 0, fmt_54, 0 }; static cilist io___51 = { 0, 3, 0, fmt_6, 0 }; static cilist io___52 = { 0, 3, 0, fmt_691, 0 }; static cilist io___54 = { 0, 3, 0, fmt_54, 0 }; static cilist io___55 = { 0, 3, 0, fmt_431, 0 }; static cilist io___58 = { 0, 6, 0, fmt_431, 0 }; static cilist io___59 = { 0, 3, 0, fmt_54, 0 }; static cilist io___60 = { 0, 6, 0, fmt_54, 0 }; static cilist io___61 = { 0, 5, 0, 0, 0 }; static cilist io___63 = { 0, 3, 0, fmt_6, 0 }; static cilist io___64 = { 0, 3, 0, fmt_66, 0 }; static cilist io___65 = { 0, 3, 0, fmt_2, 0 }; static cilist io___66 = { 0, 6, 0, fmt_2, 0 }; static cilist io___67 = { 0, 3, 0, fmt_6, 0 }; static cilist io___68 = { 0, 3, 0, fmt_21, 0 }; static cilist io___70 = { 0, 3, 0, fmt_6, 0 }; static cilist io___71 = { 0, 3, 0, fmt_6, 0 }; static cilist io___72 = { 0, 3, 0, fmt_21, 0 }; static cilist io___74 = { 0, 5, 0, 0, 0 }; static cilist io___77 = { 0, 6, 0, 0, 0 }; static cilist io___78 = { 0, 5, 0, 0, 0 }; static cilist io___81 = { 0, 6, 0, 0, 0 }; static cilist io___82 = { 0, 5, 0, 0, 0 }; static cilist io___84 = { 0, 6, 0, 0, 0 }; static cilist io___85 = { 0, 5, 0, 0, 0 }; static cilist io___87 = { 0, 6, 0, 0, 0 }; static cilist io___88 = { 0, 5, 0, 0, 0 }; static cilist io___91 = { 0, 6, 0, 0, 0 }; static cilist io___94 = { 0, 6, 0, 0, 0 }; static cilist io___96 = { 0, 6, 0, 0, 0 }; static cilist io___97 = { 0, 6, 0, 0, 0 }; static cilist io___98 = { 0, 6, 0, 0, 0 }; static cilist io___112 = { 0, 6, 0, 0, 0 }; static cilist io___113 = { 0, 6, 0, 0, 0 }; static cilist io___114 = { 0, 6, 0, 0, 0 }; static cilist io___115 = { 0, 6, 0, 0, 0 }; static cilist io___121 = { 0, 6, 0, 0, 0 }; static cilist io___122 = { 0, 6, 0, 0, 0 }; static cilist io___123 = { 0, 6, 0, fmt_43, 0 }; static cilist io___124 = { 0, 6, 0, fmt_54, 0 }; static cilist io___125 = { 0, 6, 0, 0, 0 }; static cilist io___128 = { 0, 6, 0, fmt_43, 0 }; static cilist io___130 = { 0, 6, 0, fmt_54, 0 }; static cilist io___135 = { 0, 3, 0, fmt_421, 0 }; static cilist io___136 = { 0, 3, 0, fmt_691, 0 }; static cilist io___137 = { 0, 3, 0, fmt_54, 0 }; static cilist io___138 = { 0, 3, 0, fmt_431, 0 }; static cilist io___139 = { 0, 3, 0, fmt_54, 0 }; static cilist io___140 = { 0, 3, 0, fmt_42, 0 }; static cilist io___141 = { 0, 3, 0, fmt_691, 0 }; static cilist io___142 = { 0, 3, 0, fmt_54, 0 }; static cilist io___143 = { 0, 3, 0, fmt_431, 0 }; static cilist io___144 = { 0, 3, 0, fmt_54, 0 }; static cilist io___145 = { 0, 3, 0, fmt_421, 0 }; static cilist io___146 = { 0, 3, 0, fmt_691, 0 }; static cilist io___147 = { 0, 3, 0, fmt_54, 0 }; static cilist io___148 = { 0, 3, 0, fmt_431, 0 }; static cilist io___149 = { 0, 3, 0, fmt_54, 0 }; static cilist io___150 = { 0, 3, 0, fmt_42, 0 }; static cilist io___151 = { 0, 3, 0, fmt_691, 0 }; static cilist io___152 = { 0, 3, 0, fmt_54, 0 }; static cilist io___153 = { 0, 3, 0, fmt_431, 0 }; static cilist io___154 = { 0, 3, 0, fmt_54, 0 }; static cilist io___155 = { 0, 6, 0, 0, 0 }; static cilist io___156 = { 0, 6, 0, fmt_63, 0 }; static cilist io___157 = { 0, 6, 0, fmt_54, 0 }; static cilist io___158 = { 0, 6, 0, 0, 0 }; static cilist io___159 = { 0, 6, 0, fmt_63, 0 }; static cilist io___160 = { 0, 6, 0, fmt_54, 0 }; /* This program uses inverse filtering to solve for the ARMA */ /* parameters. It requires an ARMA model for initialization. */ /* The AR model parameters are used for initialization of the */ /* inverse filtering criterion minimization to estimate */ /* the AR parameters. The MA parameters are then estimated */ /* from the AR-compensated sequence using cumulant matching. */ /* This program is currently set up to estimate at most an AR24MA8 */ /* model for a maximum of 2773 equations. For a larger model, */ /* it would be necessary to increase ns3lg and nc3lg. */ /* parameter(mp=1441) */ /* parameter(mp=4873) */ /* parameter(mp=7993) */ /* Parameter mp also in subroutine fcumb3 and mincum3 (npt) */ /* ******************************************* */ /* Determine if diagonistic on impulse response to be run */ s_rsle(&io___1); do_lio(&c__9, &c__1, diag_1.mdltst, (ftnlen)5); e_rsle(); s_wsle(&io___2); do_lio(&c__9, &c__1, "Image size(on a side)?", (ftnlen)22); e_wsle(); s_rsle(&io___3); do_lio(&c__3, &c__1, (char *)&nrow1, (ftnlen)sizeof(integer)); e_rsle(); ncol1 = nrow1; /* ***************************************************** */ /* Data read in */ s_wsle(&io___6); do_lio(&c__9, &c__1, "Is program to run an impulse response (imrsp) \ or image (image)?", (ftnlen)75); e_wsle(); s_rsle(&io___7); do_lio(&c__9, &c__1, input, (ftnlen)20); e_rsle(); if (s_cmp(input, "imrsp", (ftnlen)20, (ftnlen)5) == 0) { imprsp2_(h__, label, &nrow1, (ftnlen)30); } else { s_wsle(&io___11); do_lio(&c__9, &c__1, "B&W Image?", (ftnlen)10); e_wsle(); s_rsle(&io___12); do_lio(&c__9, &c__1, labelim, (ftnlen)30); e_rsle(); /* path='/home/teh1m/images/'//labelim */ /* call srdhead(path,string,nrow2,ncol2,iflr) */ /* Read in image */ /* do 10 j=1,ncol2 */ /* call rdcol(xr(1,j),nrow2,iflr) */ /* 10 continue */ /* call closefl(iflr) */ /* print *, 'image dimensions=',nrow2,ncol2 */ } s_wsle(&io___14); do_lio(&c__9, &c__1, "Label for predictor coefficient array?", (ftnlen)38) ; e_wsle(); s_rsle(&io___15); do_lio(&c__9, &c__1, label, (ftnlen)30); e_rsle(); s_wsle(&io___16); do_lio(&c__9, &c__1, "Number of Monte Carlo runs", (ftnlen)26); e_wsle(); s_rsle(&io___17); do_lio(&c__3, &c__1, (char *)&nmc, (ftnlen)sizeof(integer)); e_rsle(); s_wsle(&io___19); do_lio(&c__3, &c__1, (char *)&nmc, (ftnlen)sizeof(integer)); e_wsle(); s_wsle(&io___20); do_lio(&c__9, &c__1, "Unbiased or biased estimator?", (ftnlen)29); e_wsle(); s_rsle(&io___21); do_lio(&c__9, &c__1, est, (ftnlen)5); e_rsle(); s_wsle(&io___23); do_lio(&c__9, &c__1, est, (ftnlen)5); e_wsle(); s_wsle(&io___24); do_lio(&c__9, &c__1, "Table of random number generator seeds?", (ftnlen) 39); e_wsle(); s_rsle(&io___25); do_lio(&c__9, &c__1, sdtbl, (ftnlen)20); e_rsle(); s_wsle(&io___27); do_lio(&c__9, &c__1, sdtbl, (ftnlen)20); e_wsle(); /* path='/home/teh1m/suprgn/'//sdtbl */ /* Writing concatenation */ i__1[0] = 2, a__1[0] = "./"; i__1[1] = 20, a__1[1] = sdtbl; s_cat(path, a__1, i__1, &c__2, (ftnlen)80); o__1.oerr = 0; o__1.ounit = 3; o__1.ofnmlen = 80; o__1.ofnm = path; o__1.orl = 0; o__1.osta = "old"; o__1.oacc = "sequential"; o__1.ofm = "formatted"; o__1.oblnk = 0; f_open(&o__1); s_rsfe(&io___29); do_fio(&c__1, string, (ftnlen)80); e_rsfe(); s_rsfe(&io___31); do_fio(&c__1, (char *)&nsd, (ftnlen)sizeof(integer)); e_rsfe(); s_rsfe(&io___33); i__2 = nsd; for (i__ = 1; i__ <= i__2; ++i__) { do_fio(&c__1, (char *)&sd[i__ - 1], (ftnlen)sizeof(integer)); } e_rsfe(); cl__1.cerr = 0; cl__1.cunit = 3; cl__1.csta = 0; f_clos(&cl__1); /* **************************************** */ /* Inserting ARMA model parameters for initialization of minimization?' */ s_rsle(&io___36); do_lio(&c__9, &c__1, label1, (ftnlen)30); do_lio(&c__3, &c__1, (char *)&pset1, (ftnlen)sizeof(integer)); e_rsle(); s_wsle(&io___39); do_lio(&c__9, &c__1, label1, (ftnlen)30); do_lio(&c__3, &c__1, (char *)&pset1, (ftnlen)sizeof(integer)); e_wsle(); /* path='/home/teh1m/predcff/'//'cff'//label1 */ /* Writing concatenation */ i__1[0] = 5, a__1[0] = "./cff"; i__1[1] = 30, a__1[1] = label1; s_cat(path, a__1, i__1, &c__2, (ftnlen)80); o__1.oerr = 0; o__1.ounit = 3; o__1.ofnmlen = 80; o__1.ofnm = path; o__1.orl = 0; o__1.osta = 0; o__1.oacc = "sequential"; o__1.ofm = "formatted"; o__1.oblnk = 0; f_open(&o__1); al__1.aerr = 0; al__1.aunit = 3; f_rew(&al__1); s_rsfe(&io___40); do_fio(&c__1, string, (ftnlen)80); e_rsfe(); s_rsfe(&io___41); do_fio(&c__1, (char *)&dat2_1.nar, (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&nset, (ftnlen)sizeof(integer)); e_rsfe(); s_rsfe(&io___43); e_rsfe(); ncol = 1; L551: i__2 = dat2_1.nar; for (i__ = 1; i__ <= i__2; ++i__) { s_rsfe(&io___45); for (j = 1; j <= 2; ++j) { do_fio(&c__1, (char *)&dat2_1.arlst[j + (i__ << 2) - 5], (ftnlen) sizeof(integer)); } do_fio(&c__1, (char *)&ar1[i__ + ncol * 25 - 26], (ftnlen)sizeof( doublereal)); e_rsfe(); s_wsfe(&io___48); for (j = 1; j <= 2; ++j) { do_fio(&c__1, (char *)&dat2_1.arlst[j + (i__ << 2) - 5], (ftnlen) sizeof(integer)); } do_fio(&c__1, (char *)&ar1[i__ + ncol * 25 - 26], (ftnlen)sizeof( doublereal)); e_wsfe(); /* L413: */ } s_rsfe(&io___49); e_rsfe(); s_wsfe(&io___50); e_wsfe(); ++ncol; if (ncol <= nset) { goto L551; } s_rsfe(&io___51); do_fio(&c__1, string, (ftnlen)80); e_rsfe(); s_rsfe(&io___52); do_fio(&c__1, (char *)&nma0, (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&nset, (ftnlen)sizeof(integer)); e_rsfe(); s_rsfe(&io___54); e_rsfe(); ncol = 1; L553: i__2 = nma0; for (i__ = 1; i__ <= i__2; ++i__) { s_rsfe(&io___55); for (j = 1; j <= 2; ++j) { do_fio(&c__1, (char *)&malst0[j + (i__ << 2) - 5], (ftnlen)sizeof( integer)); } do_fio(&c__1, (char *)&ma0[i__ + ncol * 25 - 26], (ftnlen)sizeof( doublereal)); e_rsfe(); s_wsfe(&io___58); for (j = 1; j <= 2; ++j) { do_fio(&c__1, (char *)&malst0[j + (i__ << 2) - 5], (ftnlen)sizeof( integer)); } do_fio(&c__1, (char *)&ma0[i__ + ncol * 25 - 26], (ftnlen)sizeof( doublereal)); e_wsfe(); /* L415: */ } s_rsfe(&io___59); e_rsfe(); s_wsfe(&io___60); e_wsfe(); ++ncol; if (ncol <= nset) { goto L553; } cl__1.cerr = 0; cl__1.cunit = 3; cl__1.csta = 0; f_clos(&cl__1); /* ************************************************ */ /* Allow for the following MA support in the Inverse Filtering */ /* criterion and the stopping condition. */ s_rsle(&io___61); do_lio(&c__9, &c__1, label3, (ftnlen)30); e_rsle(); /* path='/home/teh1m/suprgn/'//label3 */ /* Writing concatenation */ i__1[0] = 2, a__1[0] = "./"; i__1[1] = 30, a__1[1] = label3; s_cat(path, a__1, i__1, &c__2, (ftnlen)80); o__1.oerr = 0; o__1.ounit = 3; o__1.ofnmlen = 80; o__1.ofnm = path; o__1.orl = 0; o__1.osta = 0; o__1.oacc = "sequential"; o__1.ofm = "formatted"; o__1.oblnk = 0; f_open(&o__1); s_rsfe(&io___63); do_fio(&c__1, string, (ftnlen)80); e_rsfe(); s_rsfe(&io___64); do_fio(&c__1, (char *)&dat2_1.nma, (ftnlen)sizeof(integer)); e_rsfe(); i__2 = dat2_1.nma; for (j = 1; j <= i__2; ++j) { s_rsfe(&io___65); for (i__ = 1; i__ <= 2; ++i__) { do_fio(&c__1, (char *)&dat2_1.malst[i__ + (j << 2) - 5], (ftnlen) sizeof(integer)); } e_rsfe(); s_wsfe(&io___66); for (i__ = 1; i__ <= 2; ++i__) { do_fio(&c__1, (char *)&dat2_1.malst[i__ + (j << 2) - 5], (ftnlen) sizeof(integer)); } e_wsfe(); /* L55: */ } s_rsfe(&io___67); do_fio(&c__1, string, (ftnlen)80); e_rsfe(); s_rsfe(&io___68); for (j = 1; j <= 2; ++j) { for (i__ = 1; i__ <= 2; ++i__) { do_fio(&c__1, (char *)&q[i__ + (j << 1) - 3], (ftnlen)sizeof( integer)); } } e_rsfe(); s_rsfe(&io___70); do_fio(&c__1, string, (ftnlen)80); e_rsfe(); s_rsfe(&io___71); do_fio(&c__1, string, (ftnlen)80); e_rsfe(); s_rsfe(&io___72); for (j = 1; j <= 2; ++j) { for (i__ = 1; i__ <= 2; ++i__) { do_fio(&c__1, (char *)&c__[i__ + (j << 1) - 3], (ftnlen)sizeof( integer)); } } e_rsfe(); cl__1.cerr = 0; cl__1.cunit = 3; cl__1.csta = 0; f_clos(&cl__1); /* ********************************************************* */ /* Resulting SNR after adding noise */ s_rsle(&io___74); do_lio(&c__3, &c__1, (char *)&ndb, (ftnlen)sizeof(integer)); do_lio(&c__9, &c__1, adnse, (ftnlen)5); e_rsle(); s_wsle(&io___77); do_lio(&c__3, &c__1, (char *)&ndb, (ftnlen)sizeof(integer)); do_lio(&c__9, &c__1, adnse, (ftnlen)5); e_wsle(); /* Blurring image flag with the following model */ s_rsle(&io___78); do_lio(&c__9, &c__1, blurflg, (ftnlen)20); do_lio(&c__9, &c__1, label3, (ftnlen)30); do_lio(&c__3, &c__1, (char *)&pset3, (ftnlen)sizeof(integer)); e_rsle(); s_wsle(&io___81); do_lio(&c__9, &c__1, blurflg, (ftnlen)20); do_lio(&c__9, &c__1, label3, (ftnlen)30); do_lio(&c__3, &c__1, (char *)&pset3, (ftnlen)sizeof(integer)); e_wsle(); /* Set flag for MA estiamtion */ s_rsle(&io___82); do_lio(&c__9, &c__1, maflg, (ftnlen)20); e_rsle(); s_wsle(&io___84); do_lio(&c__9, &c__1, maflg, (ftnlen)20); e_wsle(); /* Set flag to write parameters as file of feature vectors */ s_rsle(&io___85); do_lio(&c__9, &c__1, ftrflg, (ftnlen)10); e_rsle(); s_wsle(&io___87); do_lio(&c__9, &c__1, ftrflg, (ftnlen)10); e_wsle(); /* MA cumulant matching weighting parameter and stopping condition */ s_rsle(&io___88); do_lio(&c__5, &c__1, (char *)>ol1, (ftnlen)sizeof(doublereal)); e_rsle(); /* ********************************************************** */ /* Generate coordinates for partitioned image */ if (s_cmp(input, "image", (ftnlen)20, (ftnlen)5) == 0) { /* l=nrow2/nrow1 */ } else { l = 1; } if (l % 2 != 0) { s_wsle(&io___91); do_lio(&c__9, &c__1, "Record size is not power of sample (subimage) \ size", (ftnlen)50); e_wsle(); } k = 0; i__2 = l - 1; for (j = 0; j <= i__2; ++j) { i__3 = l - 1; for (i__ = 0; i__ <= i__3; ++i__) { ++k; cnr[(k << 1) - 2] = i__ * nrow1; cnr[(k << 1) - 1] = j * ncol1; s_wsle(&io___94); do_lio(&c__9, &c__1, "corner=", (ftnlen)7); do_lio(&c__3, &c__1, (char *)&cnr[(k << 1) - 2], (ftnlen)sizeof( integer)); do_lio(&c__3, &c__1, (char *)&cnr[(k << 1) - 1], (ftnlen)sizeof( integer)); e_wsle(); /* L28: */ } } if (s_cmp(input, "image", (ftnlen)20, (ftnlen)5) == 0) { nmc = k; } nspl = k; s_wsle(&io___96); do_lio(&c__9, &c__1, "Number of subimages=", (ftnlen)20); do_lio(&c__3, &c__1, (char *)&nspl, (ftnlen)sizeof(integer)); e_wsle(); s_wsle(&io___97); do_lio(&c__9, &c__1, "Number of initial parameter sets=", (ftnlen)33); do_lio(&c__3, &c__1, (char *)&nset, (ftnlen)sizeof(integer)); e_wsle(); if (pset1 == 0 && nmc != nset) { s_wsle(&io___98); do_lio(&c__9, &c__1, "Number of I.C. sets do not match \ number of runs", (ftnlen)80); e_wsle(); goto L680; } /* **************************************************** */ /* Begin Monte Carlo Loop */ i__3 = nmc; for (i8 = 1; i8 <= i__3; ++i8) { /* Image read in */ if (s_cmp(input, "imrsp", (ftnlen)20, (ftnlen)5) == 0) { /* Generate synthetic image */ s_copy(nse, "exp", (ftnlen)5, (ftnlen)3); iseed = sd[i8 - 1]; /* iseed=sd(25) */ imgmdl2_(h__, x, &nrow1, &iseed, nse, (ftnlen)5); /* go to 624 */ uxy = (float)0.; i__2 = ncol1; for (il2 = 1; il2 <= i__2; ++il2) { i__4 = nrow1; for (ik2 = 1; ik2 <= i__4; ++ik2) { ik3 = ik2 + 20; il3 = il2 + 20; nrf = x[ik2 + (il2 << 9) - 513]; xy[ik3 + il3 * 552 - 553] = nrf; uxy += xy[ik3 + il3 * 552 - 553]; /* L217: */ } } if (s_cmp(diag_1.mdltst, "y", (ftnlen)5, (ftnlen)1) == 0) { uxy = (float)0.; } else { uxy /= nrow1 * ncol1; } xyvar = (float)0.; i__4 = ncol1; for (il2 = 1; il2 <= i__4; ++il2) { i__2 = nrow1; for (ik2 = 1; ik2 <= i__2; ++ik2) { ik3 = ik2 + 20; il3 = il2 + 20; xy[ik3 + il3 * 552 - 553] -= uxy; /* Computing 2nd power */ d__1 = xy[ik3 + il3 * 552 - 553]; xyvar += d__1 * d__1; /* L227: */ } } if (s_cmp(diag_1.mdltst, "y", (ftnlen)5, (ftnlen)1) == 0) { xyvar = (float)1.; } else { xyvar /= nrow1 * ncol1; } } else { uxy = (float)0.; i__2 = ncol1; for (il2 = 1; il2 <= i__2; ++il2) { i__4 = nrow1; for (ik2 = 1; ik2 <= i__4; ++ik2) { xy[ik2 + 20 + (il2 + 20) * 552 - 553] = (real) xr[ik2 + cnr[(i8 << 1) - 2] + (il2 + cnr[(i8 << 1) - 1] << 9) - 513]; uxy += xy[ik2 + 20 + (il2 + 20) * 552 - 553]; /* L219: */ } } uxy /= nrow1 * ncol1; xyvar = (float)0.; i__4 = ncol1; for (il2 = 1; il2 <= i__4; ++il2) { i__2 = nrow1; for (ik2 = 1; ik2 <= i__2; ++ik2) { xy[ik2 + 20 + (il2 + 20) * 552 - 553] -= uxy; /* Computing 2nd power */ d__1 = xy[ik2 + 20 + (il2 + 20) * 552 - 553]; xyvar += d__1 * d__1; /* L220: */ } } xyvar /= nrow1 * ncol1; } s_wsle(&io___112); do_lio(&c__9, &c__1, "nrow1=", (ftnlen)6); do_lio(&c__3, &c__1, (char *)&nrow1, (ftnlen)sizeof(integer)); e_wsle(); s_wsle(&io___113); do_lio(&c__9, &c__1, "ncol1=", (ftnlen)6); do_lio(&c__3, &c__1, (char *)&ncol1, (ftnlen)sizeof(integer)); e_wsle(); s_wsle(&io___114); do_lio(&c__9, &c__1, "sample image mean=", (ftnlen)18); do_lio(&c__5, &c__1, (char *)&uxy, (ftnlen)sizeof(doublereal)); e_wsle(); s_wsle(&io___115); do_lio(&c__9, &c__1, "sample image variance=", (ftnlen)22); do_lio(&c__5, &c__1, (char *)&xyvar, (ftnlen)sizeof(doublereal)); e_wsle(); /* ************************************************* */ /* Make copy of image to generate AR-compensated sequence */ i__2 = ncol1; for (il2 = 1; il2 <= i__2; ++il2) { i__4 = nrow1; for (ik2 = 1; ik2 <= i__4; ++ik2) { xy1[ik2 + 20 + (il2 + 20) * 552 - 553] = xy[ik2 + 20 + (il2 + 20) * 552 - 553]; /* L221: */ } } /* ************************************************************** */ /* Estimation of the AR parameters by inverse filtering */ /* print *, 'initial parameters for minimization' */ /* Use a set of I.C.s or a single I.C. */ if (pset1 == 0) { pset2 = i8; } else { pset2 = pset1; } /* Initial AR estimates */ i__4 = dat2_1.nar; for (j = 1; j <= i__4; ++j) { /* if ((input.eq.'imrsp').and.(i8.eq.1)) then */ /* ar1(j,pset2)=ar1(j,pset2)+(0.10)*((-1.0)**j) */ /* end if */ p[j - 1] = ar1[j + pset2 * 25 - 26]; /* write (*,789) arlst(1,j),arlst(2,j),p(j) */ /* L791: */ } dfinv1_(p, xy, &nrow1, &ncol1, >ol1); i__4 = dat2_1.nar; for (i__ = 1; i__ <= i__4; ++i__) { ar2[i__ - 1] = p[i__ - 1]; ar[i__ + ((i8 << 2) + 1) * 25 - 126] = ar2[i__ - 1]; /* write (*,789) arlst(1,i),arlst(2,i),ar2(i) */ /* L789: */ /* L93: */ } /* **************************************************** */ /* Generate AR compensated sequence */ /* do 222 l1=1,ncol1 */ /* do 222 k1=1,nrow1 */ /* l2=l1+20 */ /* k2=k1+20 */ /* buf1=0.0 */ /* do 224 k=1,nar */ /* i2=arlst(1,k) */ /* j2=arlst(2,k) */ /* buf1=buf1+ar2(k)*xy1(k2-i2,l2-j2) */ /* 224 continue */ /* xy(k2,l2)=buf1 */ /* 222 continue */ /* Check center of impulse response array */ /* ictr=nrow1/2+20 */ /* print *, 'Impulse response' */ /* write (*,533) ((xy1(i,j),j=ictr-2,ictr+4),i=ictr-2,ictr+4) */ /* print *, 'AR compensated sequence' */ /* write (*,533) ((xy(i,j),j=ictr-2,ictr+4),i=ictr-2,ictr+4) */ /* 533 format (7e11.3) */ /* ******************************************* */ /* Estimate MA parameters using the C(q,k-c) algorithm */ /* call cqk(p,xy,nrow1,ncol1,q,c) */ /* print *, 'MA estimates from C(q,k-c) alg.' */ /* do 96 i=1,nma */ /* write (*,789) malst(1,i),malst(2,i),p(i) */ /* 96 continue */ /* stop */ /* ******************************************* */ /* Estimate gamma2 and gamma3 */ /* do 732 j=1,nma */ /* ma1(j)=p(j) */ /* 732 continue */ /* call gamma_MA(ma1,xy,nrow1,ncol1,est,gm2,gm3) */ /* ************************************************************ */ /* if (MAflg.eq.'MAest') then */ /* go to 751 */ /* if (mdltst.eq.'y') then */ /* do 734 j=1,nma-1 */ /* p(j)=p(j)+ ((-1.0)**j)*(0.2) */ /* 734 continue */ /* end if */ /* do 94 i=1,nma */ /* ma(i,1,i8)=p(i) */ /* 94 continue */ /* End of MA estimation */ /* ****************************************** */ s_wsle(&io___121); do_lio(&c__9, &c__1, "Monte Carlo run number", (ftnlen)22); do_lio(&c__3, &c__1, (char *)&i8, (ftnlen)sizeof(integer)); e_wsle(); s_wsle(&io___122); do_lio(&c__9, &c__1, "The AR parameters are", (ftnlen)21); e_wsle(); i__4 = dat2_1.nar; for (i__ = 1; i__ <= i__4; ++i__) { s_wsfe(&io___123); for (j = 1; j <= 2; ++j) { do_fio(&c__1, (char *)&dat2_1.arlst[j + (i__ << 2) - 5], ( ftnlen)sizeof(integer)); } do_fio(&c__1, (char *)&ar[i__ + ((i8 << 2) + 1) * 25 - 126], ( ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&ar1[i__ + pset2 * 25 - 26], (ftnlen)sizeof( doublereal)); e_wsfe(); /* L419: */ } s_wsfe(&io___124); e_wsfe(); s_wsle(&io___125); do_lio(&c__9, &c__1, "The MA parameters are", (ftnlen)21); e_wsle(); if (s_cmp(ftrgam, "gamma2", (ftnlen)30, (ftnlen)6) == 0) { /* ma(nma,1,i8)=gm2/xyvar */ } else { ma[dat2_1.nma + ((i8 << 2) + 1) * 25 - 126] = (float)1.; } i__4 = dat2_1.nma; for (i__ = 1; i__ <= i__4; ++i__) { s_wsfe(&io___128); for (j = 1; j <= 2; ++j) { do_fio(&c__1, (char *)&dat2_1.malst[j + (i__ << 2) - 5], ( ftnlen)sizeof(integer)); } do_fio(&c__1, (char *)&ma[i__ + ((i8 << 2) + 1) * 25 - 126], ( ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&ma1[i__ - 1], (ftnlen)sizeof(doublereal)); e_wsfe(); /* L418: */ } s_wsfe(&io___130); e_wsfe(); /* L624: */ /* End Monte Carlo loop */ /* L625: */ } /* *************************************************************** */ /* Calculate statistics of the Monte Carlo Runs */ for (i__ = 1; i__ <= 1; ++i__) { i__3 = dat2_1.nar; for (j = 1; j <= i__3; ++j) { i__4 = nmc; for (i8 = 1; i8 <= i__4; ++i8) { /* L720: */ uar[j + i__ * 25 - 26] += ar[j + (i__ + (i8 << 2)) * 25 - 126] ; } } } for (i__ = 1; i__ <= 1; ++i__) { i__4 = dat2_1.nma; for (j = 1; j <= i__4; ++j) { i__3 = nmc; for (i8 = 1; i8 <= i__3; ++i8) { /* L722: */ uma[j + i__ * 25 - 26] += ma[j + (i__ + (i8 << 2)) * 25 - 126] ; } } } for (i__ = 1; i__ <= 1; ++i__) { i__3 = dat2_1.nar; for (j = 1; j <= i__3; ++j) { uar[j + i__ * 25 - 26] /= nmc; sdar[j + i__ * 25 - 26] = (float)0.; i__4 = nmc; for (i8 = 1; i8 <= i__4; ++i8) { /* L726: */ /* Computing 2nd power */ d__1 = ar[j + (i__ + (i8 << 2)) * 25 - 126] - uar[j + i__ * 25 - 26]; sdar[j + i__ * 25 - 26] += d__1 * d__1; } sdar[j + i__ * 25 - 26] = sqrt(sdar[j + i__ * 25 - 26] / nmc); /* L725: */ } } for (i__ = 1; i__ <= 1; ++i__) { i__3 = dat2_1.nma; for (j = 1; j <= i__3; ++j) { uma[j + i__ * 25 - 26] /= nmc; sdma[j + i__ * 25 - 26] = (float)0.; i__4 = nmc; for (i8 = 1; i8 <= i__4; ++i8) { /* L728: */ /* Computing 2nd power */ d__1 = ma[j + (i__ + (i8 << 2)) * 25 - 126] - uma[j + i__ * 25 - 26]; sdma[j + i__ * 25 - 26] += d__1 * d__1; } sdma[j + i__ * 25 - 26] = sqrt(sdma[j + i__ * 25 - 26] / nmc); /* L727: */ } } /* ********************************************************** */ if (s_cmp(ftrflg, "ftron", (ftnlen)10, (ftnlen)5) != 0) { nmc = 1; /* Store Mean ARMA model coefficients */ /* path='/home/teh1m/predcff/'//'cff'//label */ /* Writing concatenation */ i__1[0] = 5, a__1[0] = "./cff"; i__1[1] = 30, a__1[1] = label; s_cat(path, a__1, i__1, &c__2, (ftnlen)80); o__1.oerr = 0; o__1.ounit = 3; o__1.ofnmlen = 80; o__1.ofnm = path; o__1.orl = 0; o__1.osta = 0; o__1.oacc = "sequential"; o__1.ofm = "formatted"; o__1.oblnk = 0; f_open(&o__1); al__1.aerr = 0; al__1.aunit = 3; f_rew(&al__1); s_wsfe(&io___135); e_wsfe(); s_wsfe(&io___136); do_fio(&c__1, (char *)&dat2_1.nar, (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&nmc, (ftnlen)sizeof(integer)); e_wsfe(); s_wsfe(&io___137); e_wsfe(); ncol = 1; L548: i__3 = dat2_1.nar; for (i__ = 1; i__ <= i__3; ++i__) { s_wsfe(&io___138); for (j = 1; j <= 2; ++j) { do_fio(&c__1, (char *)&dat2_1.arlst[j + (i__ << 2) - 5], ( ftnlen)sizeof(integer)); } do_fio(&c__1, (char *)&uar[i__ + ncol * 25 - 26], (ftnlen)sizeof( doublereal)); e_wsfe(); /* L412: */ } s_wsfe(&io___139); e_wsfe(); ++ncol; if (ncol <= nmc) { goto L548; } s_wsfe(&io___140); e_wsfe(); s_wsfe(&io___141); do_fio(&c__1, (char *)&dat2_1.nma, (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&nmc, (ftnlen)sizeof(integer)); e_wsfe(); s_wsfe(&io___142); e_wsfe(); ncol = 1; L549: i__3 = dat2_1.nma; for (i__ = 1; i__ <= i__3; ++i__) { s_wsfe(&io___143); for (j = 1; j <= 2; ++j) { do_fio(&c__1, (char *)&dat2_1.malst[j + (i__ << 2) - 5], ( ftnlen)sizeof(integer)); } do_fio(&c__1, (char *)&uma[i__ + ncol * 25 - 26], (ftnlen)sizeof( doublereal)); e_wsfe(); /* L41: */ } s_wsfe(&io___144); e_wsfe(); ++ncol; if (ncol <= nmc) { goto L549; } cl__1.cerr = 0; cl__1.cunit = 3; cl__1.csta = 0; f_clos(&cl__1); } /* ************************************* */ /* Store parameters as feature vectors */ if (s_cmp(ftrflg, "ftron", (ftnlen)10, (ftnlen)5) == 0) { /* Store parameter feature vectors */ /* path='/home/teh1m/predcff/cff'//label */ /* Writing concatenation */ i__1[0] = 5, a__1[0] = "./cff"; i__1[1] = 30, a__1[1] = label; s_cat(path, a__1, i__1, &c__2, (ftnlen)80); o__1.oerr = 0; o__1.ounit = 3; o__1.ofnmlen = 80; o__1.ofnm = path; o__1.orl = 0; o__1.osta = 0; o__1.oacc = "sequential"; o__1.ofm = "formatted"; o__1.oblnk = 0; f_open(&o__1); al__1.aerr = 0; al__1.aunit = 3; f_rew(&al__1); s_wsfe(&io___145); e_wsfe(); s_wsfe(&io___146); do_fio(&c__1, (char *)&dat2_1.nar, (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&nmc, (ftnlen)sizeof(integer)); e_wsfe(); s_wsfe(&io___147); e_wsfe(); ncol = 1; L648: i__3 = dat2_1.nar; for (i__ = 1; i__ <= i__3; ++i__) { s_wsfe(&io___148); for (j = 1; j <= 2; ++j) { do_fio(&c__1, (char *)&dat2_1.arlst[j + (i__ << 2) - 5], ( ftnlen)sizeof(integer)); } do_fio(&c__1, (char *)&ar[i__ + ((ncol << 2) + 1) * 25 - 126], ( ftnlen)sizeof(doublereal)); e_wsfe(); /* L612: */ } s_wsfe(&io___149); e_wsfe(); ++ncol; if (ncol <= nmc) { goto L648; } s_wsfe(&io___150); e_wsfe(); s_wsfe(&io___151); do_fio(&c__1, (char *)&dat2_1.nma, (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&nmc, (ftnlen)sizeof(integer)); e_wsfe(); s_wsfe(&io___152); e_wsfe(); ncol = 1; L649: i__3 = dat2_1.nma; for (i__ = 1; i__ <= i__3; ++i__) { s_wsfe(&io___153); for (j = 1; j <= 2; ++j) { do_fio(&c__1, (char *)&dat2_1.malst[j + (i__ << 2) - 5], ( ftnlen)sizeof(integer)); } do_fio(&c__1, (char *)&ma[i__ + ((ncol << 2) + 1) * 25 - 126], ( ftnlen)sizeof(doublereal)); e_wsfe(); /* L31: */ } s_wsfe(&io___154); e_wsfe(); ++ncol; if (ncol <= nmc) { goto L649; } /* L523: */ cl__1.cerr = 0; cl__1.cunit = 3; cl__1.csta = 0; f_clos(&cl__1); } /* **************************************************** */ /* Display the mean paramester estimates */ s_wsle(&io___155); do_lio(&c__9, &c__1, "The AR mean parameters are", (ftnlen)26); e_wsle(); /* L545: */ i__3 = dat2_1.nar; for (i__ = 1; i__ <= i__3; ++i__) { s_wsfe(&io___156); for (j = 1; j <= 2; ++j) { do_fio(&c__1, (char *)&dat2_1.arlst[j + (i__ << 2) - 5], (ftnlen) sizeof(integer)); } do_fio(&c__1, (char *)&uar[i__ - 1], (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&sdar[i__ - 1], (ftnlen)sizeof(doublereal)); e_wsfe(); /* L541: */ } s_wsfe(&io___157); e_wsfe(); s_wsle(&io___158); do_lio(&c__9, &c__1, "The MA mean parameters are", (ftnlen)26); e_wsle(); /* L547: */ i__3 = dat2_1.nma; for (i__ = 1; i__ <= i__3; ++i__) { s_wsfe(&io___159); for (j = 1; j <= 2; ++j) { do_fio(&c__1, (char *)&dat2_1.malst[j + (i__ << 2) - 5], (ftnlen) sizeof(integer)); } do_fio(&c__1, (char *)&uma[i__ - 1], (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&sdma[i__ - 1], (ftnlen)sizeof(doublereal)); e_wsfe(); /* L542: */ } s_wsfe(&io___160); e_wsfe(); L680: ; } /* MAIN__ */ /* Subroutine */ int gamma_ma__(ma1, xy, nrow1, ncol1, est, gm2, gm3, est_len) doublereal *ma1, *xy; integer *nrow1, *ncol1; char *est; doublereal *gm2, *gm3; ftnlen est_len; { /* System generated locals */ integer i__1, i__2; doublereal d__1; /* Builtin functions */ integer s_cmp(), s_wsle(), do_lio(), e_wsle(); /* Local variables */ static integer i__, j; static doublereal b20, b30, buf2, buf3; /* Fortran I/O blocks */ static cilist io___167 = { 0, 6, 0, 0, 0 }; static cilist io___168 = { 0, 6, 0, 0, 0 }; /* Parameter adjustments */ xy -= 553; --ma1; /* Function Body */ b30 = (float)0.; b20 = (float)0.; i__1 = dat2_1.nma; for (j = 1; j <= i__1; ++j) { /* Computing 3rd power */ d__1 = ma1[j]; b30 += d__1 * (d__1 * d__1); /* Computing 2nd power */ d__1 = ma1[j]; b20 += d__1 * d__1; /* L82: */ } buf3 = (float)0.; buf2 = (float)0.; i__1 = *ncol1 + 20; for (j = 21; j <= i__1; ++j) { i__2 = *nrow1 + 20; for (i__ = 21; i__ <= i__2; ++i__) { buf3 += xy[i__ + j * 552] * xy[i__ + j * 552] * xy[i__ + j * 552]; buf2 += xy[i__ + j * 552] * xy[i__ + j * 552]; /* L23: */ } } if (s_cmp(diag_1.mdltst, "y", (ftnlen)5, (ftnlen)1) != 0) { buf2 /= *nrow1 * *ncol1; buf3 /= *nrow1 * *ncol1; } *gm3 = buf3 / b30; *gm2 = buf2 / b20; s_wsle(&io___167); do_lio(&c__9, &c__1, "Estimate of gamma2=", (ftnlen)19); do_lio(&c__5, &c__1, (char *)&(*gm2), (ftnlen)sizeof(doublereal)); e_wsle(); s_wsle(&io___168); do_lio(&c__9, &c__1, "Estimate of gamma3=", (ftnlen)19); do_lio(&c__5, &c__1, (char *)&(*gm3), (ftnlen)sizeof(doublereal)); e_wsle(); /* L91: */ return 0; } /* gamma_ma__ */ /* Main program alias */ int armaest22_ () { MAIN__ (); }