diff options
Diffstat (limited to 'noao/imred/dtoi/hdicfit')
31 files changed, 3453 insertions, 0 deletions
diff --git a/noao/imred/dtoi/hdicfit/hdic.com b/noao/imred/dtoi/hdicfit/hdic.com new file mode 100644 index 00000000..959df10e --- /dev/null +++ b/noao/imred/dtoi/hdicfit/hdic.com @@ -0,0 +1,6 @@ +# Common block for hdtoi package/ icgfit package interface. + +int nraw +double maxden +pointer den, big_den +common /raw/maxden, den, big_den, nraw diff --git a/noao/imred/dtoi/hdicfit/hdicadd.x b/noao/imred/dtoi/hdicfit/hdicadd.x new file mode 100644 index 00000000..1fae152d --- /dev/null +++ b/noao/imred/dtoi/hdicfit/hdicadd.x @@ -0,0 +1,47 @@ +include "hdicfit.h" + +# HDIC_ADDPOINT -- Add a density, exposure, weights point into the sample. + +procedure hdic_addpoint (ic, + nden, nexp, nwts, den, y, wts, userwts, x, wd, sdev, npts) + +pointer ic # Pointer to ic data structure +real nden # New density value to be added +real nexp # New exposure value to be added +double nwts # New weight value to be added +pointer den # Pointer to existing density array +pointer y # Pointer to existing exposure array +pointer wts # Pointer to existing wts array +pointer userwts # Pointer to existing userwts array +pointer x # Pointer to existing array of ind variables +pointer wd # Pointer to flag array for deletion reasons +pointer sdev # Pointer to standard deviation array +int npts # Number of points, incremented on output + +begin + npts = npts + 1 + + call realloc (den, npts, TY_DOUBLE) + call realloc (y, npts, TY_DOUBLE) + call realloc (wts, npts, TY_DOUBLE) + call realloc (userwts, npts, TY_DOUBLE) + call realloc (x, npts, TY_DOUBLE) + call realloc (wd, npts, TY_INT) + call realloc (sdev, npts, TY_DOUBLE) + + Memd[den+npts-1] = double (nden) + Memd[y +npts-1] = double (nexp) + Memd[wts+npts-1] = double (nwts) + Memd[userwts+npts-1] = double (nwts) + Memi[wd +npts-1] = NDELETE + Memd[sdev+npts-1] = ADDED_PT + + # Sort the data and then update the reference vector. + call hdic_sort (Memd[den], Memd[y], Memd[wts], Memd[userwts], + Memi[wd], Memd[sdev], npts) + call hdic_init (Memd[den], npts, Memd[den+npts-1]) + + IC_NEWX(ic) = YES + IC_NEWY(ic) = YES + IC_NEWWTS(ic) = YES +end diff --git a/noao/imred/dtoi/hdicfit/hdicclean.x b/noao/imred/dtoi/hdicfit/hdicclean.x new file mode 100644 index 00000000..6c838123 --- /dev/null +++ b/noao/imred/dtoi/hdicfit/hdicclean.x @@ -0,0 +1,94 @@ +include <pkg/rg.h> +include "hdicfit.h" + +# IC_CLEAN -- Replace rejected points by the fitted values. + +procedure ic_cleand (ic, cv, x, y, w, npts) + +pointer ic # ICFIT pointer +pointer cv # Curfit pointer +double x[ARB] # Ordinates +double y[ARB] # Abscissas +double w[ARB] # Weights +int npts # Number of points + +int i, nclean, newreject +pointer sp, xclean, yclean, wclean +double dcveval() + +begin + if ((IC_LOW(ic) == 0.) && (IC_HIGH(ic) == 0.)) + return + + # If there has been no subsampling and no sample averaging then the + # IC_REJPTS(ic) array already contains the rejected points. + + if (npts == IC_NFIT(ic)) { + if (IC_NREJECT(ic) > 0) { + do i = 1, npts { + if (Memi[IC_REJPTS(ic)+i-1] == YES) + y[i] = dcveval (cv, x[i]) + } + } + + # If there has been no sample averaging then the rejpts array already + # contains indices into the subsampled array. + + } else if (abs(IC_NAVERAGE(ic)) == 1) { + if (IC_NREJECT(ic) > 0) { + do i = 1, npts { + if (Memi[IC_REJPTS(ic)+i-1] == YES) + Memd[IC_YFIT(ic)+i-1] = + dcveval (cv, Memd[IC_XFIT(ic)+i-1]) + } + } + call rg_unpackd (IC_RG(ic), Memd[IC_YFIT(ic)], y) + + # Because ic_fit rejects points from the fitting data which + # has been sample averaged the rejpts array refers to the wrong data. + # Do the cleaning using ic_deviant to find the points to reject. + + } else if (RG_NPTS(IC_RG(ic)) == npts) { + call amovki (NO, Memi[IC_REJPTS(ic)], npts) + call ic_deviantd (cv, x, y, w, Memi[IC_REJPTS(ic)], npts, + IC_LOW(ic), IC_HIGH(ic), IC_GROW(ic), NO, IC_NREJECT(ic), + newreject) + if (IC_NREJECT(ic) > 0) { + do i = 1, npts { + if (Memi[IC_REJPTS(ic)+i-1] == YES) + y[i] = dcveval (cv, x[i]) + } + } + + # If there is subsampling then allocate temporary arrays for the + # subsample points. + + } else { + call smark (sp) + + nclean = RG_NPTS(IC_RG(ic)) + call salloc (xclean, nclean, TY_DOUBLE) + call salloc (yclean, nclean, TY_DOUBLE) + call salloc (wclean, nclean, TY_DOUBLE) + + call rg_packd (IC_RG(ic), x, Memd[xclean]) + call rg_packd (IC_RG(ic), y, Memd[yclean]) + call rg_packd (IC_RG(ic), w, Memd[wclean]) + call amovki (NO, Memi[IC_REJPTS(ic)], npts) + + call ic_deviantd (cv, Memd[xclean], Memd[yclean], + Memd[wclean], Memi[IC_REJPTS(ic)], nclean, IC_LOW(ic), + IC_HIGH(ic), IC_GROW(ic), NO, IC_NREJECT(ic), newreject) + + if (IC_NREJECT(ic) > 0) { + do i = 1, npts { + if (Memi[IC_REJPTS(ic)+i-1] == YES) + Memd[yclean+i-1] = dcveval (cv, Memd[xclean+i-1]) + } + } + call rg_unpackd (IC_RG(ic), Memd[yclean], y) + + call sfree (sp) + } +end + diff --git a/noao/imred/dtoi/hdicfit/hdicdeviant.x b/noao/imred/dtoi/hdicfit/hdicdeviant.x new file mode 100644 index 00000000..0c6f8f75 --- /dev/null +++ b/noao/imred/dtoi/hdicfit/hdicdeviant.x @@ -0,0 +1,116 @@ +include <mach.h> +include <math/curfit.h> + +# IC_DEVIANT -- Find deviant points with large residuals from the fit +# and reject from the fit. +# +# The sigma of the fit residuals is calculated. The rejection thresholds +# are set at +-reject*sigma. Points outside the rejection threshold are +# recorded in the reject array. + +procedure ic_deviantd (cv, x, y, w, rejpts, npts, low_reject, high_reject, + grow, refit, nreject, newreject) + +pointer cv # Curve descriptor +double x[ARB] # Input ordinates +double y[ARB] # Input data values +double w[ARB] # Weights +int rejpts[ARB] # Points rejected +int npts # Number of input points +real low_reject, high_reject # Rejection thresholds +real grow # Rejection radius +int refit # Refit the curve? +int nreject # Number of points rejected +int newreject # Number of new points rejected + +int i, j, i_min, i_max, ilast +double sigma, low_cut, high_cut, residual +pointer sp, residuals + +begin + # If low_reject and high_reject are zero then simply return. + if ((low_reject == 0.) && (high_reject == 0.)) + return + + # Allocate memory for the residuals. + call smark (sp) + call salloc (residuals, npts, TY_DOUBLE) + + # Compute the residuals. + call dcvvector (cv, x, Memd[residuals], npts) + call asubd (y, Memd[residuals], Memd[residuals], npts) + + # Compute the sigma of the residuals. If there are less than + # 5 points return. + + j = 0 + nreject = 0 + sigma = 0. + + do i = 1, npts { + if ((w[i] != 0.) && (rejpts[i] == NO)) { + sigma = sigma + Memd[residuals+i-1] ** 2 + j = j + 1 + } else if (rejpts[i] == YES) + nreject = nreject + 1 + } + + if (j < 5) { + call sfree (sp) + return + } else + sigma = sqrt (sigma / j) + + if (low_reject > 0.) + low_cut = -low_reject * sigma + else + low_cut = -MAX_REAL + if (high_reject > 0.) + high_cut = high_reject * sigma + else + high_cut = MAX_REAL + + # Reject the residuals exceeding the rejection limits. + # A for loop is used instead of do because with region growing we + # want to modify the loop index. + + newreject = 0 + for (i = 1; i <= npts; i = i + 1) { + if ((w[i] == 0.) || (rejpts[i] == YES)) + next + + residual = Memd[residuals + i - 1] + if ((residual > high_cut) || (residual < low_cut)) { + i_min = max (1, int (i - grow)) + i_max = min (npts, int (i + grow)) + + # Reject points from the fit and flag them. + do j = i_min, i_max { + if ((abs (x[i] - x[j]) <= grow) && (w[j] != 0.) && + (rejpts[j] == NO)) { + if (refit == YES) + call dcvrject (cv, x[j], y[j], w[j]) + rejpts[j] = YES + newreject = newreject + 1 + ilast = j + } + } + i = ilast + } + } + + if ((refit == YES) && (newreject > 0)) { + call dcvsolve (cv, i) + + switch (i) { + case SINGULAR: + call eprintf ("ic_reject: Singular solution\n") + case NO_DEG_FREEDOM: + call eprintf ("ic_reject: No degrees of freedom\n") + } + } + + nreject = nreject + newreject + + call sfree (sp) +end diff --git a/noao/imred/dtoi/hdicfit/hdicdosetup.x b/noao/imred/dtoi/hdicfit/hdicdosetup.x new file mode 100644 index 00000000..c9e117ec --- /dev/null +++ b/noao/imred/dtoi/hdicfit/hdicdosetup.x @@ -0,0 +1,104 @@ +include <math/curfit.h> +include "hdicfit.h" + +# IC_DOSETUP -- Setup the fit. This is called at the start of each call +# to ic_fit to update the fitting parameters if necessary. + +procedure ic_dosetupd (ic, cv, x, wts, npts, newx, newwts, newfunction, refit) + +pointer ic # ICFIT pointer +pointer cv # Curfit pointer +double x[ARB] # Ordinates of data +double wts[ARB] # Weights +int npts # Number of points in data +int newx # New x points? +int newwts # New weights? +int newfunction # New function? +int refit # Use cvrefit? + +int ord +pointer rg_xrangesd() +extern hd_powerd() +errchk rg_xrangesd, malloc + +begin + # Set sample points. + if ((newx == YES) || (newwts == YES)) { + if (npts == 0) + call error (0, "No data points for fit") + + call mfree (IC_XFIT(ic), TY_DOUBLE) + call mfree (IC_YFIT(ic), TY_DOUBLE) + call malloc (IC_XFIT(ic), npts, TY_DOUBLE) + + call mfree (IC_WTSFIT(ic), TY_DOUBLE) + call malloc (IC_WTSFIT(ic), npts, TY_DOUBLE) + + call mfree (IC_REJPTS(ic), TY_INT) + call malloc (IC_REJPTS(ic), npts, TY_INT) + call amovki (NO, Memi[IC_REJPTS(ic)], npts) + IC_NREJECT(ic) = 0 + + # Set sample points. + + call rg_free (IC_RG(ic)) + IC_RG(ic) = rg_xrangesd (Memc[IC_SAMPLE(ic)], x, npts) + call rg_order (IC_RG(ic)) + call rg_merge (IC_RG(ic)) + call rg_wtbind (IC_RG(ic), abs (IC_NAVERAGE(ic)), x, wts, npts, + Memd[IC_XFIT(ic)], Memd[IC_WTSFIT(ic)], IC_NFIT(ic)) + + if (IC_NFIT(ic) == 0) + call error (0, "No sample points for fit") + + if (IC_NFIT(ic) == npts) { + call rg_free (IC_RG(ic)) + call mfree (IC_XFIT(ic), TY_DOUBLE) + call mfree (IC_WTSFIT(ic), TY_DOUBLE) + IC_YFIT(ic) = NULL + IC_WTSFIT(ic) = NULL + } else + call malloc (IC_YFIT(ic), IC_NFIT(ic), TY_DOUBLE) + + refit = NO + } + + # Set curve fitting parameters. + + if ((newx == YES) || (newfunction == YES)) { + if (cv != NULL) + call dcvfree (cv) + + switch (IC_FUNCTION(ic)) { + case LEGENDRE, CHEBYSHEV: + ord = min (IC_ORDER(ic), IC_NFIT(ic)) + call dcvinit (cv, IC_FUNCTION(ic), ord, double (IC_XMIN(ic)), + double (IC_XMAX(ic))) + case SPLINE1: + ord = min (IC_ORDER(ic), IC_NFIT(ic) - 1) + if (ord > 0) + call dcvinit (cv, SPLINE1, ord, double (IC_XMIN(ic)), + double (IC_XMAX(ic))) + else + call dcvinit (cv, LEGENDRE, IC_NFIT(ic), + double (IC_XMIN(ic)), double (IC_XMAX(ic))) + case SPLINE3: + ord = min (IC_ORDER(ic), IC_NFIT(ic) - 3) + if (ord > 0) + call dcvinit (cv, SPLINE3, ord, double (IC_XMIN(ic)), + double (IC_XMAX(ic))) + else + call dcvinit (cv, LEGENDRE, IC_NFIT(ic), + double (IC_XMIN(ic)), double (IC_XMAX(ic))) + case USERFNC: + ord = min (IC_ORDER(ic), IC_NFIT(ic)) + call dcvinit (cv, USERFNC, ord, double (IC_XMIN(ic)), + double (IC_XMAX(ic))) + call dcvuserfnc (cv, hd_powerd) + default: + call error (0, "Unknown fitting function") + } + + refit = NO + } +end diff --git a/noao/imred/dtoi/hdicfit/hdicebars.x b/noao/imred/dtoi/hdicfit/hdicebars.x new file mode 100644 index 00000000..335c161d --- /dev/null +++ b/noao/imred/dtoi/hdicfit/hdicebars.x @@ -0,0 +1,217 @@ +include <gset.h> +include <mach.h> +include <pkg/gtools.h> +include "hdicfit.h" + +define MSIZE 2.0 + +# HDIC_EBARS -- Plot data points with their error bars as markers. The +# independent variable is fog subtracted transformed density. + +procedure hdic_ebars (ic, gp, gt, x, y, wts, ebw, npts) + +pointer ic # Pointer to ic structure +pointer gp # Pointer to graphics stream +pointer gt # Pointer to gtools structure +double x[ARB] # Array of independent variables +double y[ARB] # Array of dependent variables +double wts[ARB] # Array of weights +double ebw[ARB] # Error bar half width (positive number) +int npts # Number of points + +pointer sp, xr, yr, sz +int orig_mark, i, xaxis, yaxis, mark, added_mark +real size, dx, dy, szmk +int gt_geti() +bool fp_equald() +include "hdic.com" + +begin + call smark (sp) + + xaxis = (IC_AXES (ic, IC_GKEY(ic), 1)) + yaxis = (IC_AXES (ic, IC_GKEY(ic), 2)) + + if (xaxis == 'x' || xaxis == 'u') + orig_mark = GM_HEBAR + else if (yaxis == 'x' || yaxis == 'u') + orig_mark = GM_VEBAR + else { + call eprintf ("Choose graph type with axes 'u' or 'x'\n") + call sfree (sp) + return + } + + added_mark = GM_CIRCLE + + call salloc (xr, npts, TY_REAL) + call salloc (yr, npts, TY_REAL) + call salloc (sz, npts, TY_REAL) + + call achtdr (x, Memr[xr], npts) + call achtdr (y, Memr[yr], npts) + call achtdr (ebw, Memr[sz], npts) + + if (IC_OVERPLOT(ic) == NO) { + # Start a new plot + call gclear (gp) + + # Set the graph scale and axes + call gascale (gp, Memr[xr], npts, 1) + call gascale (gp, Memr[yr], npts, 2) + + # If plotting HD curve, set wy2 to maxden, which may have + # been updated if a new endpoint was added. + + if ((IC_AXES (ic, IC_GKEY(ic), 1) == 'y') && + (IC_AXES (ic, IC_GKEY(ic), 2) == 'u')) + call gswind (gp, INDEF, INDEF, INDEF, real (maxden)) + + call gt_swind (gp, gt) + call gt_labax (gp, gt) + } + + call ggscale (gp, 0.0, 0.0, dx, dy) + + do i = 1, npts { + size = Memr[sz+i-1] # Sizes are WCS units; transform them + if (gt_geti (gt, GTTRANSPOSE) == NO) { + size = size / dx + mark = orig_mark + szmk = size + # Check for added point + if (fp_equald (ebw[i], ADDED_PT)) { + szmk = MSIZE + mark = added_mark + } + + # Check for deleted point + if (fp_equald (wts[i], 0.0D0)) { + szmk = MSIZE + mark = mark + GM_CROSS + } + + call gmark (gp, Memr[xr+i-1], Memr[yr+i-1], mark, szmk, szmk) + + } else { + size = size / dy + szmk = size + mark = orig_mark + + # Check for added point + if (fp_equald (ebw[i], ADDED_PT)) { + szmk = MSIZE + mark = added_mark + } + + # Check for deleted point + if (fp_equald (wts[i], 0.0D0)) { + szmk = MSIZE + mark = mark + GM_CROSS + } + + call gmark (gp, Memr[yr+i-1], Memr[xr+i-1], mark, szmk, szmk) + } + } + + IC_OVERPLOT(ic) = NO + call sfree (sp) +end + + +# HDIC_EBW -- Calculate error bar width for plotting points. Width is +# returned in NDC units. + +procedure hdic_ebw (ic, density, indv, sdev, ebw, npts) + +pointer ic # Pointer to ic structure +double density[ARB] # Untransformed density NOT fog subtracted +double indv[ARB] # Transformed density above fog +double sdev[ARB] # Array of standard deviation values, density units +double ebw[ARB] # Error bar half width (positive numbers) +int npts # Number of data points + +double fog +pointer sp, denaf, outwe +int xaxis, yaxis, i +bool fp_equald() +real ic_getr() + +begin + xaxis = (IC_AXES (ic, IC_GKEY(ic), 1)) + yaxis = (IC_AXES (ic, IC_GKEY(ic), 2)) + + if (xaxis == 'u' || yaxis == 'u') { + call amovd (sdev, ebw, npts) + return + } + + call smark (sp) + call salloc (denaf, npts, TY_DOUBLE) + call salloc (outwe, npts, TY_DOUBLE) + + fog = double (ic_getr (ic, "fog")) + + call asubkd (density, fog, Memd[denaf], npts) + call aaddd (Memd[denaf], sdev, Memd[denaf], npts) + + call hdic_ebtran (Memd[denaf], Memd[outwe], npts, IC_TRANSFORM(ic)) + + # Subtract transformed values to get errors. Then check for + # added points, which are flagged with ebw=ADDED_PT. + + call asubd (Memd[outwe], indv, ebw, npts) + do i = 1, npts { + if (fp_equald (sdev[i], ADDED_PT)) + ebw[i] = ADDED_PT + } + + call sfree (sp) +end + + +# HDIC_EBTRAN -- Apply transformation, generating a vector of independent +# variables from a density vector. The independent variable vector is the +# standard deviation values and the output values are the errors. This +# routine checks for ADDED_PT valued standard deviation, which indicates an +# added point. + +procedure hdic_ebtran (density, ind_var, npts, transform) + +double density[npts] # Density vector - input +double ind_var[npts] # Ind variable vector - filled on output +int npts # Length of data vectors +int transform # Integer code for transform type + +int i +bool fp_equald() + +begin + switch (transform) { + case HD_LOGO: + do i = 1, npts { + if (fp_equald (density[i], 0.0)) + ind_var[i] = 0.0 + else + ind_var[i] = log10 ((10. ** density[i]) - 1.0) + } + case HD_K75: + do i = 1, npts { + if (fp_equald (density[i], 0.0)) + ind_var[i] = 0.0 + else + ind_var[i] = density[i] + 0.75*log10(1.- 10.**(-density[i])) + } + case HD_K50: + do i = 1, npts { + if (fp_equald (density[i], 0.0)) + ind_var[i] = 0.0 + else + ind_var[i] = density[i] + 0.50*log10(1.- 10.**(-density[i])) + } + case HD_NONE: + call amovd (density, ind_var, npts) + default: + call error (0, "Unrecognized transformation in HDIC_EBTRAN") + } +end diff --git a/noao/imred/dtoi/hdicfit/hdicerrors.x b/noao/imred/dtoi/hdicfit/hdicerrors.x new file mode 100644 index 00000000..7722e5e2 --- /dev/null +++ b/noao/imred/dtoi/hdicfit/hdicerrors.x @@ -0,0 +1,143 @@ +include <math/curfit.h> +include "hdicfit.h" + +# IC_ERRORS -- Compute and error diagnositic information. + +procedure ic_errorsd (ic, file, cv, x, y, wts, npts) + +pointer ic # ICFIT pointer +char file[ARB] # Output file +pointer cv # Curfit pointer +double x[ARB] # Ordinates +double y[ARB] # Abscissas +double wts[ARB] # Weights +int npts # Number of data points + +int i, n, deleted, ncoeffs, fd +double chisqr, rms +pointer sp, fit, wts1, coeffs, errors + +int dcvstati(), open() +double ic_rmsd() +errchk open() + +begin + # Open the output file. + fd = open (file, APPEND, TEXT_FILE) + + # Determine the number of coefficients and allocate memory. + ncoeffs = dcvstati (cv, CVNCOEFF) + call smark (sp) + call salloc (coeffs, ncoeffs, TY_DOUBLE) + call salloc (errors, ncoeffs, TY_DOUBLE) + + if (npts == IC_NFIT(ic)) { + # Allocate memory for the fit. + n = npts + call salloc (fit, n, TY_DOUBLE) + call salloc (wts1, n, TY_DOUBLE) + + # Eliminate rejected points and count deleted points. + call amovd (wts, Memd[wts1], n) + if (IC_NREJECT(ic) > 0) { + do i = 1, npts { + if (Memi[IC_REJPTS(ic)+i-1] == YES) + Memd[wts1+i-1] = 0. + } + } + deleted = 0 + do i = 1, n { + if (wts[i] == 0.) + deleted = deleted + 1 + } + + # Get the coefficients and compute the errors. + call dcvvector (cv, x, Memd[fit], n) + call dcvcoeff (cv, Memd[coeffs], ncoeffs) + call dcverrors (cv, y, Memd[wts1], Memd[fit], n, chisqr, + Memd[errors]) + rms = ic_rmsd (x, y, Memd[fit], Memd[wts1], n) + + } else { + # Allocate memory for the fit. + n = IC_NFIT(ic) + call salloc (fit, n, TY_DOUBLE) + call salloc (wts1, n, TY_DOUBLE) + + # Eliminate rejected points and count deleted points. + call amovd (Memd[IC_WTSFIT(ic)], Memd[wts1], n) + if (IC_NREJECT(ic) > 0) { + do i = 1, npts { + if (Memi[IC_REJPTS(ic)+i-1] == YES) + Memd[wts1+i-1] = 0. + } + } + deleted = 0 + do i = 1, n { + if (wts[i] == 0.) + deleted = deleted + 1 + } + + # Get the coefficients and compute the errors. + call dcvvector (cv, Memd[IC_XFIT(ic)], Memd[fit], n) + rms = ic_rmsd (Memd[IC_XFIT(ic)], Memd[IC_YFIT(ic)], + Memd[fit], Memd[wts1], n) + call dcvcoeff (cv, Memd[coeffs], ncoeffs) + call dcverrors (cv, Memd[IC_YFIT(ic)], Memd[wts1], Memd[fit], + n, chisqr, Memd[errors]) + } + + # Print the error analysis. + call fprintf (fd, "total points = %d\nsample points = %d\n") + call pargi (npts) + call pargi (n) + call fprintf (fd, "nrejected = %d\ndeleted = %d\n") + call pargi (IC_NREJECT(ic)) + call pargi (deleted) + call fprintf (fd, "RMS = %7.4g\n") + call pargd (rms) + call fprintf (fd, "square root of reduced chi square = %7.4g\n") + call pargd (sqrt (chisqr)) + + #call fprintf (fd, "\tcoefficent\terror\n") + #do i = 1, ncoeffs { + # call fprintf (fd, "\t%10.4e\t%10.4e\n") + # call parg$t (Mem$t[coeffs+i-1]) + # call parg$t (Mem$t[errors+i-1]) + #} + + # Free allocated memory. + call sfree (sp) + call close (fd) +end + + +# IC_RMS -- Compute RMS of points which have not been deleted. + +double procedure ic_rmsd (x, y, fit, wts, npts) + +double x[ARB] # Ordinates +double y[ARB] # Abscissas +double fit[ARB] # Fit +double wts[ARB] # Weights +int npts # Number of data points + +int i, n +double resid, rms + +begin + rms = 0. + n = 0 + do i = 1, npts { + if (wts[i] == 0.) + next + resid = y[i] - fit[i] + rms = rms + resid * resid + n = n + 1 + } + + if (n > 0) + rms = sqrt (rms / n) + + return (rms) +end diff --git a/noao/imred/dtoi/hdicfit/hdicfit.h b/noao/imred/dtoi/hdicfit/hdicfit.h new file mode 100644 index 00000000..89d9c73d --- /dev/null +++ b/noao/imred/dtoi/hdicfit/hdicfit.h @@ -0,0 +1,65 @@ +# Definition file for DTOI task hdfit which uses the hdicfit subdirectory. + +define NSPOTS 64 # Initially, max number of calibration spots +define NVALS_FIT 200 # Number of points in vector of fitted points +define WT_NONE 0 # No weighting used in fit +define WT_USER 1 # User specifies weighting in fit +define WT_CALC 2 # Weighting in fit calculated from std dev +define MIN_DEN EPSILONR # Density used for setting curfit minval +define HD_NONE 1 # Ind var is density - no transformation +define HD_LOGO 2 # Ind var is log opacitance +define HD_K50 3 # Ind var is Kaiser transform w/ alpha=0.50 +define HD_K75 4 # Ind var is Kaiser transform w/ alpha=0.75 +define UDELETE 100 # Point deleted by user flag +define PDELETE 101 # Point deleted by program +define NDELETE 102 # Point not deleted +define ADDED_PT 0.0 # Indication of added point in sdev array + +# The ICFIT data structure - modified for use with the DTOI package. + +define IC_NGKEYS 5 # Number of graph keys +define IC_LENSTRUCT 47 # Length of ICFIT structure + +# User fitting parameters +define IC_FUNCTION Memi[$1] # Function type +define IC_ORDER Memi[$1+1] # Order of function +define IC_SAMPLE Memi[$1+2] # Pointer to sample string +define IC_NAVERAGE Memi[$1+3] # Sampling averaging bin +define IC_NITERATE Memi[$1+4] # Number of rejection interation +define IC_TRANSFORM Memi[$1+5] # Type of transformation ** DTOI ONLY ** +define IC_XMIN Memr[P2R($1+6)] # Minimum value for curve +define IC_XMAX Memr[P2R($1+7)] # Maximum value for curve +define IC_LOW Memr[P2R($1+8)] # Low rejection value +define IC_HIGH Memr[P2R($1+9)] # Low rejection value +define IC_GROW Memr[P2R($1+10)]# Rejection growing radius + +# ICFIT parameters used for fitting +define IC_NFIT Memi[$1+11] # Number of fit points +define IC_NREJECT Memi[$1+12] # Number of rejected points +define IC_RG Memi[$1+13] # Pointer for ranges +define IC_XFIT Memi[$1+14] # Pointer to ordinates of fit points +define IC_YFIT Memi[$1+15] # Pointer to abscissas of fit points +define IC_WTSFIT Memi[$1+16] # Pointer to weights of fit points +define IC_REJPTS Memi[$1+17] # Pointer to rejected points + +# ICFIT parameters used for interactive graphics +define IC_NEWX Memi[$1+18] # New x fit points? +define IC_NEWY Memi[$1+19] # New y points? +define IC_NEWWTS Memi[$1+20] # New weights? +define IC_NEWFUNCTION Memi[$1+21] # New fitting function? +define IC_NEWTRANSFORM Memi[$1+22] # New transform? ** DTOI ONLY ** +define IC_OVERPLOT Memi[$1+23] # Overplot next plot? +define IC_FITERROR Memi[$1+24] # Error in fit +define IC_LABELS Memi[$1+25+$2-1]# Graph axis labels +define IC_UNITS Memi[$1+27+$2-1]# Graph axis units + +define IC_FOG Memr[P2R($1+29)]# *** DTOI ONLY *** value of fog level +define IC_NEWFOG Memi[$1+30] # Flag for change in fog +define IC_RESET Memi[$1+31] # Flag for resetting variables +define IC_UPDATE Memi[$1+32] +define IC_EBARS Memi[$1+33] # Flag for plotting error bars +define IC_RFOG Memr[P2R($1+34)]# Reference value of fog for resetting + +# ICFIT key definitions +define IC_GKEY Memi[$1+35] # Graph key +define IC_AXES Memi[$1+36+($2-1)*2+$3-1] # Graph axis codes diff --git a/noao/imred/dtoi/hdicfit/hdicfit.x b/noao/imred/dtoi/hdicfit/hdicfit.x new file mode 100644 index 00000000..91e05f9d --- /dev/null +++ b/noao/imred/dtoi/hdicfit/hdicfit.x @@ -0,0 +1,80 @@ +include <math/curfit.h> +include "hdicfit.h" + +# IC_FIT -- Fit a function. This is the main fitting task. It uses +# flags to define changes since the last fit. This allows the most +# efficient use of the curfit and ranges packages. + +procedure ic_fitd (ic, cv, x, y, wts, npts, newx, newy, newwts, newfunction) + +pointer ic # ICFIT pointer +pointer cv # Curfit pointer +double x[npts] # Ordinates +double y[npts] # Data to be fit +double wts[npts] # Weights +int npts # Number of points +int newx # New x points? +int newy # New y points? +int newwts # New weights? +int newfunction # New function? + +int ier, refit +errchk ic_dosetupd + +begin + # Setup the new parameters. + call ic_dosetupd (ic, cv, x, wts, npts, newx, newwts, newfunction, + refit) + + if (npts == IC_NFIT(ic)) { + # If not sampling use the data array directly. + if (refit == NO) { + call dcvfit (cv, x, y, wts, npts, WTS_USER, ier) + } else if (newy == YES) + call dcvrefit (cv, x, y, wts, ier) + + } else { + # If sampling first form the sample y values. + if ((newx == YES) || (newy == YES) || (newwts == YES)) + call rg_wtbind (IC_RG(ic), IC_NAVERAGE(ic), y, wts, npts, + Memd[IC_YFIT(ic)], Memd[IC_WTSFIT(ic)], IC_NFIT(ic)) + if (refit == NO) { + call dcvfit (cv, Memd[IC_XFIT(ic)], Memd[IC_YFIT(ic)], + Memd[IC_WTSFIT(ic)], IC_NFIT(ic), WTS_USER, ier) + } else if (newy == YES) + call dcvrefit (cv, Memd[IC_XFIT(ic)], Memd[IC_YFIT(ic)], + Memd[IC_WTSFIT(ic)], ier) + } + + # Check for an error in the fit. + switch (ier) { + case SINGULAR: + call printf ("Singular solution") + call flush (STDOUT) + case NO_DEG_FREEDOM: + call printf ("No degrees of freedom") + call flush (STDOUT) + return + } + + refit = YES + + # Do pixel rejection if desired. + if ((IC_LOW(ic) > 0.) || (IC_HIGH(ic) > 0.)) { + if (npts == IC_NFIT(ic)) { + call ic_rejectd (cv, x, y, wts, Memi[IC_REJPTS(ic)], + IC_NFIT(ic), IC_LOW(ic), IC_HIGH(ic), IC_NITERATE(ic), + IC_GROW(ic), IC_NREJECT(ic)) + } else { + call ic_rejectd (cv, Memd[IC_XFIT(ic)], Memd[IC_YFIT(ic)], + Memd[IC_WTSFIT(ic)], Memi[IC_REJPTS(ic)], IC_NFIT(ic), + IC_LOW(ic), IC_HIGH(ic), IC_NITERATE(ic), IC_GROW(ic), + IC_NREJECT(ic)) + } + + if (IC_NREJECT(ic) > 0) + refit = NO + } else + IC_NREJECT(ic) = 0 +end + diff --git a/noao/imred/dtoi/hdicfit/hdicgaxes.x b/noao/imred/dtoi/hdicfit/hdicgaxes.x new file mode 100644 index 00000000..4f016c9d --- /dev/null +++ b/noao/imred/dtoi/hdicfit/hdicgaxes.x @@ -0,0 +1,101 @@ +include <pkg/gtools.h> +include "hdicfit.h" + +# ICG_AXES -- Set axes data. +# The applications program may set additional axes types. + +procedure icg_axesd (ic, gt, cv, axis, x, y, z, npts) + +pointer ic # ICFIT pointer +pointer gt # GTOOLS pointer +pointer cv # CURFIT pointer +int axis # Output axis +double x[npts] # Independent variable +double y[npts] # Dependent variable +double z[npts] # Output values +int npts # Number of points + +int i, axistype, gtlabel[2], gtunits[2] +double a, b, xmin, xmax +pointer label, units + +double dcveval(), icg_dvzd() +errchk adivd() +extern icg_dvzd() + +data gtlabel/GTXLABEL, GTYLABEL/ +data gtunits/GTXUNITS, GTYUNITS/ + +begin + axistype = IC_AXES(ic, IC_GKEY(ic), axis) + switch (axistype) { + case 'x': # Independent variable + call gt_sets (gt, gtlabel[axis], Memc[IC_LABELS(ic,1)]) + call gt_sets (gt, gtunits[axis], Memc[IC_UNITS(ic,1)]) + call amovd (x, z, npts) + call gt_sets (gt, GTSUBTITLE, "") + case 'y': # Dependent variable + call gt_sets (gt, gtlabel[axis], Memc[IC_LABELS(ic,2)]) + call gt_sets (gt, gtunits[axis], Memc[IC_UNITS(ic,2)]) + call amovd (y, z, npts) + call gt_sets (gt, GTSUBTITLE, "") + case 'f': # Fitted values + call gt_sets (gt, gtlabel[axis], "fit") + call gt_sets (gt, gtunits[axis], Memc[IC_UNITS(ic,2)]) + call dcvvector (cv, x, z, npts) + call gt_sets (gt, GTSUBTITLE, "") + case 'r': # Residuals + call gt_sets (gt, gtlabel[axis], "residuals") + call gt_sets (gt, gtunits[axis], Memc[IC_UNITS(ic,2)]) + call dcvvector (cv, x, z, npts) + call asubd (y, z, z, npts) + call gt_sets (gt, GTSUBTITLE, "") + case 'd': # Ratio + call gt_sets (gt, gtlabel[axis], "ratio") + call gt_sets (gt, gtunits[axis], "") + call dcvvector (cv, x, z, npts) +# iferr (call adiv$t (y, z, z, npts)) + call advzd (y, z, z, npts, icg_dvzd) + call gt_sets (gt, GTSUBTITLE, "") + case 'n': # Linear component removed + call gt_sets (gt, gtlabel[axis], "non-linear component") + call gt_sets (gt, gtunits[axis], Memc[IC_UNITS(ic,2)]) + xmin = IC_XMIN(ic) + xmax = IC_XMAX(ic) + a = dcveval (cv, double(xmin)) + b = (dcveval (cv, double(xmax)) - a) / (xmax - xmin) + do i = 1, npts + z[i] = y[i] - a - b * (x[i] - xmin) + call gt_sets (gt, GTSUBTITLE, "") + default: # User axes types. + call malloc (label, SZ_LINE, TY_CHAR) + call malloc (units, SZ_LINE, TY_CHAR) + if (axis == 1) { + call strcpy (Memc[IC_LABELS(ic,1)], Memc[label], SZ_LINE) + call strcpy (Memc[IC_UNITS(ic,1)], Memc[units], SZ_LINE) + call amovd (x, z, npts) + } else { + call strcpy (Memc[IC_LABELS(ic,2)], Memc[label], SZ_LINE) + call strcpy (Memc[IC_UNITS(ic,2)], Memc[units], SZ_LINE) + call amovd (y, z, npts) + } + call icg_uaxesd (ic, axistype, cv, x, y, z, npts, Memc[label], + Memc[units], SZ_LINE) + call gt_sets (gt, gtlabel[axis], Memc[label]) + call gt_sets (gt, gtunits[axis], Memc[units]) + call gt_sets (gt, GTSUBTITLE, "HD Curve") + call mfree (label, TY_CHAR) + call mfree (units, TY_CHAR) + } +end + + +# ICG_DVZ -- Error action to take on zero division. + +double procedure icg_dvzd (x) + +double x # Numerator + +begin + return (1.) +end diff --git a/noao/imred/dtoi/hdicfit/hdicgcolon.x b/noao/imred/dtoi/hdicfit/hdicgcolon.x new file mode 100644 index 00000000..52299df7 --- /dev/null +++ b/noao/imred/dtoi/hdicfit/hdicgcolon.x @@ -0,0 +1,284 @@ +include <error.h> +include <gset.h> +include "hdicfit.h" + +define EB_WTS 10 +define EB_SDEV 11 + +# List of colon commands. +define CMDS "|show|sample|naverage|function|order|low_reject|high_reject|\ + |niterate|grow|errors|vshow|transform|fog|reset|quit|ebars|" + +define SHOW 1 # Show values of parameters +define SAMPLE 2 # Set or show sample ranges +define NAVERAGE 3 # Set or show sample averaging or medianing +define FUNCTION 4 # Set or show function type +define ORDER 5 # Set or show order +define LOW_REJECT 6 # Set or show lower rejection factor +define HIGH_REJECT 7 # Set or show upper rejection factor +# newline 8 +define NITERATE 9 # Set or show rejection iterations +define GROW 10 # Set or show rejection growing radius +define ERRORS 11 # Show errors of fit +define VSHOW 12 # Show verbose information +define TRANSFORM 13 # Set or show transformation +define FOG 14 # Set or show value of fog +define RESET 15 # Reset x, y, wts, npts to original values +define QUIT 16 # Terminate without updating database +define EBARS 17 # Set error bars to represent weights or + # standard deviations + +# ICG_COLON -- Processes colon commands. + +procedure icg_colond (ic, cmdstr, gp, gt, cv, x, y, wts, npts) + +pointer ic # ICFIT pointer +char cmdstr[ARB] # Command string +pointer gp # GIO pointer +pointer gt # GTOOLS pointer +pointer cv # CURFIT pointer for error listing +double x[npts], y[npts], wts[npts] # Data arrays for error listing +int npts # Number of data points + +real rval +char cmd[SZ_LINE] +int ncmd, ival, ip, junk + +int nscan(), strdic(), strncmp(), ctor() +string funcs "|chebyshev|legendre|spline1|spline3|power|" +string tform "|none|logopacitance|k50|k75|" + +begin + # Use formated scan to parse the command string. + # The first word is the command and it may be minimum match + # abbreviated with the list of commands. + + call sscan (cmdstr) + call gargwrd (cmd, SZ_LINE) + ncmd = strdic (cmd, cmd, SZ_LINE, CMDS) + + switch (ncmd) { + case SHOW: + # show: Show the values of the fitting parameters. The terminal + # is cleared and paged using the gtools paging procedures. + + call gargwrd (cmd, SZ_LINE) + if (nscan() == 1) { + call gdeactivate (gp, AW_CLEAR) + call ic_show (ic, "STDOUT", gt) + call greactivate (gp, AW_PAUSE) + } else { + iferr (call ic_show (ic, cmd, gt)) + call erract (EA_WARN) + } + + case SAMPLE: + # sample: List or set the sample points. + + call gargwrd (cmd, SZ_LINE) + if (cmd[1] == EOS) { + call printf ("sample = %s\n") + call pargstr (Memc[IC_SAMPLE(ic)]) + } else { + call strcpy (cmd, Memc[IC_SAMPLE(ic)], SZ_LINE) + IC_NEWX(ic) = YES + } + + case NAVERAGE: + # naverage: List or set the sample averging. + + call gargi (ival) + if (nscan() == 1) { + call printf ("naverage = %d\n") + call pargi (IC_NAVERAGE(ic)) + } else { + IC_NAVERAGE(ic) = ival + IC_NEWX(ic) = YES + } + + case FUNCTION: + # function: List or set the fitting function. + + call gargwrd (cmd, SZ_LINE) + if (cmd[1] == EOS) { + call printf ("function = %s\n") + call ic_gstr (ic, "function", cmd, SZ_LINE) + call pargstr (cmd) + } else { + if (strdic (cmd, cmd, SZ_LINE, funcs) > 0) { + call ic_pstr (ic, "function", cmd) + IC_NEWFUNCTION(ic) = YES + } else + call printf ("Unknown or ambiguous function") + } + + case ORDER: + # order: List or set the function order. + + call gargi (ival) + if (nscan() == 1) { + call printf ("order = %d\n") + call pargi (IC_ORDER(ic)) + } else { + IC_ORDER(ic) = ival + IC_NEWFUNCTION(ic) = YES + } + + case LOW_REJECT: + # low_reject: List or set the lower rejection threshold. + + call gargr (rval) + if (nscan() == 1) { + call printf ("low_reject = %g\n") + call pargr (IC_LOW(ic)) + } else + IC_LOW(ic) = rval + + case HIGH_REJECT: + # high_reject: List or set the high rejection threshold. + + call gargr (rval) + if (nscan() == 1) { + call printf ("high_reject = %g\n") + call pargr (IC_HIGH(ic)) + } else + IC_HIGH(ic) = rval + + case NITERATE: + # niterate: List or set the number of rejection iterations. + + call gargi (ival) + if (nscan() == 1) { + call printf ("niterate = %d\n") + call pargi (IC_NITERATE(ic)) + } else + IC_NITERATE(ic) = ival + + case GROW: + # grow: List or set the rejection growing. + + call gargr (rval) + if (nscan() == 1) { + call printf ("grow = %g\n") + call pargr (IC_GROW(ic)) + } else + IC_GROW(ic) = rval + + case ERRORS: + call gargwrd (cmd, SZ_LINE) + if (nscan() == 1) { + call gdeactivate (gp, AW_CLEAR) + call ic_show (ic, "STDOUT", gt) + call ic_errorsd (ic, "STDOUT", cv, x, y, wts, npts) + call greactivate (gp, AW_PAUSE) + } else { + iferr { + call ic_show (ic, cmd, gt) + call ic_errorsd (ic, cmd, cv, x, y, wts, npts) + } then + call erract (EA_WARN) + } + case VSHOW: + # verbose show: Show the values of the fitting parameters. + # The terminal is paged using the gtools paging procedure. + + call gargwrd (cmd, SZ_LINE) + if (cmd[1] == EOS) { + call gdeactivate (gp, AW_CLEAR) + call ic_vshowd (ic, "STDOUT", cv, x, y, wts, npts, gt) + call greactivate (gp, AW_PAUSE) + } else { + iferr { + call ic_vshowd (ic, cmd, cv, x, y, wts, npts, gt) + } then + call erract (EA_WARN) + } + case TRANSFORM: + # transform: List or set the transformation type. This + # option applies to HDTOI procedures only. + + call gargwrd (cmd, SZ_LINE) + if (cmd[1] == EOS) { + call printf ("transform = %s\n") + call ic_gstr (ic, "transform", cmd, SZ_LINE) + call pargstr (cmd) + } else { + ival= strdic (cmd, cmd, SZ_LINE, tform) + if (ival > 0) { + call ic_pstr (ic, "transform", cmd) + IC_NEWTRANSFORM(ic) = YES + IC_NEWX(ic) = YES + switch (IC_TRANSFORM(ic)) { + case HD_NONE: + call ic_pstr (ic, "xlabel", "Density") + case HD_LOGO: + call ic_pstr (ic, "xlabel", + "Log Opacitance: log (10**Den - 1)") + case HD_K50: + call ic_pstr (ic, "xlabel", + "Den + 0.50 * Log (1 - (10 ** -Den))") + case HD_K75: + call ic_pstr (ic, "xlabel", + "Den + 0.75 * Log (1 - (10 ** -Den))") + } + } else + call printf ("Unknown or ambiguous transform") + } + + case FOG: + # fog: DTOI ONLY - change or reset the value of the fog level + + call gargwrd (cmd, SZ_LINE) + if (cmd[1] == EOS) { + call printf ("fog = %g\n") + call pargr (IC_FOG(ic)) + } else { + if (strncmp (cmd, "reset", 1) == 0) + IC_FOG(ic) = IC_RFOG(ic) + else { + ip = 1 + junk = ctor (cmd, ip, rval) + IC_FOG(ic) = rval + } + IC_NEWFOG(ic) = YES + IC_NEWX(ic) = YES + } + + case RESET: + # Set flag to reset x, y, wts and npts to original values. + IC_RESET(ic) = YES + IC_NEWX(ic) = YES + IC_NEWY(ic) = YES + IC_NEWWTS(ic) = YES + IC_NEWFUNCTION(ic) = YES + IC_NEWTRANSFORM(ic) = YES + + case QUIT: + # Set update flag to know + IC_UPDATE(ic) = NO + + case EBARS: + # [HV]BAR marker can indicate either errors or weights + call gargwrd (cmd, SZ_LINE) + if (cmd[1] == EOS) { + if (IC_EBARS(ic) == EB_WTS) + call printf ("ebars = Weights\n") + else if (IC_EBARS(ic) == EB_SDEV) + call printf ("ebars = Errors\n") + } else { + if (strncmp (cmd, "weights", 1) == 0 || + strncmp (cmd, "WEIGHTS", 1) == 0) + IC_EBARS(ic) = EB_WTS + else if (strncmp (cmd, "errors", 1) == 0 || + strncmp (cmd, "ERRORS", 1) == 0) + IC_EBARS(ic) = EB_SDEV + else + call printf ("Unrecognized value for ebars '%s'\n") + call pargstr (cmd) + } + + default: + call eprintf ("Unrecognized command '%s'\n") + call pargstr (cmd) + } +end diff --git a/noao/imred/dtoi/hdicfit/hdicgdelete.x b/noao/imred/dtoi/hdicfit/hdicgdelete.x new file mode 100644 index 00000000..82c61f70 --- /dev/null +++ b/noao/imred/dtoi/hdicfit/hdicgdelete.x @@ -0,0 +1,81 @@ +include <gset.h> +include <mach.h> +include <pkg/gtools.h> +include "hdicfit.h" + +define MSIZE 2. # Mark size + +# ICG_DELETE -- Delete data point nearest the cursor. +# The nearest point to the cursor in NDC coordinates is determined. + +procedure icg_deleted (ic, gp, gt, cv, x, y, wts, npts, wx, wy) + +pointer ic # ICFIT pointer +pointer gp # GIO pointer +pointer gt # GTOOLS pointer +pointer cv # CURFIT pointer +double x[npts], y[npts] # Data points +double wts[npts] # Weight array +int npts # Number of points +real wx, wy # Position to be nearest + +int gt_geti() +pointer sp, xout, yout + +begin + call smark (sp) + call salloc (xout, npts, TY_DOUBLE) + call salloc (yout, npts, TY_DOUBLE) + + call icg_axesd (ic, gt, cv, 1, x, y, Memd[xout], npts) + call icg_axesd (ic, gt, cv, 2, x, y, Memd[yout], npts) + if (gt_geti (gt, GTTRANSPOSE) == NO) + call icg_d1d (ic, gp, Memd[xout], Memd[yout], wts, npts, wx, wy) + else + call icg_d1d (ic, gp, Memd[yout], Memd[xout], wts, npts, wy, wx) + + call sfree (sp) +end + + +# ICG_D1D - Do the actual delete. + +procedure icg_d1d (ic, gp, x, y, wts, npts, wx, wy) + +pointer ic # ICFIT pointer +pointer gp # GIO pointer +double x[npts], y[npts] # Data points +double wts[npts] # Weight array +int npts # Number of points +real wx, wy # Position to be nearest + +int i, j +real x0, y0, r2, r2min + +begin + # Transform world cursor coordinates to NDC. + call gctran (gp, wx, wy, wx, wy, 1, 0) + + # Search for nearest point to a point with non-zero weight. + r2min = MAX_REAL + do i = 1, npts { + if (wts[i] == 0.) + next + + call gctran (gp, real (x[i]), real (y[i]), x0, y0, 1, 0) + + r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2 + if (r2 < r2min) { + r2min = r2 + j = i + } + } + + # Mark the deleted point with a cross and set the weight to zero. + if (j != 0) { + call gscur (gp, real (x[j]), real (y[j])) + call gmark (gp, real (x[j]), real (y[j]), GM_CROSS, MSIZE, MSIZE) + wts[j] = 0. + IC_NEWWTS(ic) = YES + } +end diff --git a/noao/imred/dtoi/hdicfit/hdicgfit.x b/noao/imred/dtoi/hdicfit/hdicgfit.x new file mode 100644 index 00000000..b50ed03c --- /dev/null +++ b/noao/imred/dtoi/hdicfit/hdicgfit.x @@ -0,0 +1,402 @@ +include <error.h> +include <pkg/gtools.h> +include <mach.h> +include <gset.h> +include "hdicfit.h" + +define HELP "noao$lib/scr/hdicgfit.key" +define PROMPT "hdicfit options" +define EB_WTS 10 +define EB_SDEV 11 + +# ICG_FIT -- Interactive curve fitting with graphics. This is the main +# entry point for the interactive graphics part of the icfit package. + +procedure icg_fitd (ic, gp, cursor, gt, cv, density, loge, owts, osdev, npts) + +pointer ic # ICFIT pointer +pointer gp # GIO pointer +char cursor[ARB] # GIO cursor input +pointer gt # GTOOLS pointer +pointer cv # CURFIT pointer +double density[npts] # Original density, not fog subtracted +double loge[npts] # Original log exposure values +double owts[npts] # Original weights array +double osdev[npts] # Original standard deviation array +int npts # Number of points + +real wx, wy, wwy, oldx, oldy +int wcs, key, nptso +char cmd[SZ_LINE] + +int i, newgraph, axes[2], linetype +double x1, newwt +pointer userwts, x, y, wts, den, whydel, sdev, ebw + +int clgcur(), stridxs(), scan(), nscan() +int icg_nearestd(), gstati() +double dcveval() +errchk ic_fitd, malloc +define redraw_ 91 + +begin + # Allocate memory for the fit and a copy of the weights. + # The weights are copied because they are changed when points are + # deleted. The input x is the untransformed density, and is used + # to generate other types of transforms. Points can be added to + # the sample, so the y array and weights array can change as well. + # The original number of points is also remembered. + + call malloc (userwts, npts, TY_DOUBLE) + call malloc (x, npts, TY_DOUBLE) + call malloc (y, npts, TY_DOUBLE) + call malloc (wts, npts, TY_DOUBLE) + call malloc (den, npts, TY_DOUBLE) + call malloc (whydel, npts, TY_INT) + call malloc (sdev, npts, TY_DOUBLE) + call malloc (ebw, npts, TY_DOUBLE) + + call amovd (owts, Memd[userwts], npts) + call amovd (owts, Memd[wts], npts) + call amovd (loge, Memd[y], npts) + call amovd (density, Memd[den], npts) + call amovki (NDELETE, Memi[whydel], npts) + call amovd (osdev, Memd[sdev], npts) + nptso = npts + + # Initialize + IC_OVERPLOT(ic) = NO + IC_NEWX(ic) = YES + IC_NEWY(ic) = YES + IC_NEWWTS(ic) = YES + IC_NEWFUNCTION(ic) = YES + IC_NEWTRANSFORM(ic) = YES + IC_UPDATE(ic) = YES + IC_EBARS(ic) = EB_SDEV + + # Read cursor commands. + + key = 'f' + axes[1] = IC_AXES(ic, IC_GKEY(ic), 1) + axes[2] = IC_AXES(ic, IC_GKEY(ic), 2) + + repeat { + switch (key) { + case '?': # Print help text. + call gpagefile (gp, HELP, PROMPT) + + case 'q': # Terminate cursor loop + break + + case ':': # List or set parameters + if (stridxs ("/", cmd) == 1) + call gt_colon (cmd, gp, gt, newgraph) + else + call icg_colond (ic, cmd, gp, gt, cv, Memd[x], + Memd[y], Memd[wts], npts) + + if (IC_RESET(ic) == YES) { + npts = nptso + call amovd (owts, Memd[userwts], npts) + call amovd (owts, Memd[wts], npts) + call amovd (loge, Memd[y], npts) + call amovd (density, Memd[den], npts) + call amovki (NDELETE, Memi[whydel], npts) + call amovd (osdev, Memd[sdev], npts) + call hdic_init (density, npts, 0.0) + } + + # See if user wants to quit without updating + if (IC_UPDATE(ic) == NO) { + call mfree (x, TY_DOUBLE) + call mfree (y, TY_DOUBLE) + call mfree (wts, TY_DOUBLE) + call mfree (userwts, TY_DOUBLE) + call mfree (den, TY_DOUBLE) + call mfree (sdev, TY_DOUBLE) + return + } + + case 'a': # Add data points to the sample. This is only possible + # from an HD curve plot. + + if ((IC_AXES (ic, IC_GKEY(ic), 1) == 'y') && + (IC_AXES (ic, IC_GKEY(ic), 2) == 'u')) { + + # Query for weight after plotting current location + # call gt_plot (gp, gt, wx, wy, 1) + call gmark (gp, wx, wy, GM_CIRCLE, 2.0, 2.0) + newwt = 1.0D0 + call printf ("Enter weight of new point (%g): ") + call pargd (newwt) + call flush (STDOUT) + if (scan() != EOF) { + call pargd (x1) + if (nscan() == 1) { + if (!IS_INDEFD (x1)) { + newwt = x1 + } + } + } + + } else { + call eprintf ("Points can be added only from an HD Curve\n") + next + } + + # Add fog into "density above fog" value read from cursor + wwy = wy + IC_FOG (ic) + if (wwy < 0.0) { + call eprintf ( + "New density (%g) is below fog and will not be added\n") + call pargr (wwy) + next + } + + # Add point into sample + call eprintf ("New Point: density above fog = %.4f, log ") + call pargr (wwy) + call eprintf ("exposure = %.4f, weight = %.4f\n") + call pargr (wx) + call pargd (newwt) + + call hdic_addpoint (ic, wwy, wx, newwt, den, y, wts, userwts, + x, whydel, sdev, npts) + + call realloc (ebw, npts, TY_DOUBLE) + + call hdic_transform (ic, Memd[den], Memd[userwts], Memd[x], + Memd[wts], Memi[whydel], npts) + + case 'c': # Print the positions of data points. + i = icg_nearestd (ic, gp, gt, cv, Memd[x], Memd[y], npts, + wx, wy) + + if (i != 0) { + call printf ("den= %7.4g x= %7.4g exp= %7.4g fit= %7.4g") + call pargd (Memd[den+i-1]) + call pargd (Memd[x+i-1]) + call pargd (Memd[y+i-1]) + call pargd (dcveval (cv, Memd[x+i-1])) + } + + case 'd': # Delete data points. + call icg_deleted (ic, gp, gt, cv, Memd[x], Memd[y], Memd[wts], + npts, wx, wy) + + case 'f': # Fit the function and reset the flags. + iferr { + # Copy new transformed vector, if necessary + if (IC_NEWTRANSFORM(ic) == YES || IC_NEWFOG(ic) == YES) + call hdic_transform (ic, Memd[den], Memd[userwts], + Memd[x], Memd[wts], Memi[whydel], npts) + + call ic_fitd (ic, cv, Memd[x], Memd[y], Memd[wts], npts, + IC_NEWX(ic), IC_NEWY(ic), IC_NEWWTS(ic), + IC_NEWFUNCTION(ic)) + + IC_NEWX(ic) = NO + IC_NEWY(ic) = NO + IC_NEWWTS(ic) = NO + IC_NEWFUNCTION(ic) = NO + IC_NEWTRANSFORM(ic) = NO + IC_FITERROR(ic) = NO + IC_NEWFOG(ic) = NO + newgraph = YES + } then { + IC_FITERROR(ic) = YES + call erract (EA_WARN) + newgraph = NO + } + + case 'g': # Set graph axes types. + call printf ("Graph key to be defined: ") + call flush (STDOUT) + if (scan() == EOF) + goto redraw_ + call gargc (cmd[1]) + + switch (cmd[1]) { + case '\n': + case 'h', 'i', 'j', 'k', 'l': + switch (cmd[1]) { + case 'h': + key = 1 + case 'i': + key = 2 + case 'j': + key = 3 + case 'k': + key = 4 + case 'l': + key = 5 + } + + call printf ("Set graph axes types (%c, %c): ") + call pargc (IC_AXES(ic, key, 1)) + call pargc (IC_AXES(ic, key, 2)) + call flush (STDOUT) + if (scan() == EOF) + goto redraw_ + call gargc (cmd[1]) + + switch (cmd[1]) { + case '\n': + default: + call gargc (cmd[2]) + call gargc (cmd[2]) + if (cmd[2] != '\n') { + IC_AXES(ic, key, 1) = cmd[1] + IC_AXES(ic, key, 2) = cmd[2] + } + } + default: + call printf ("Not a graph key") + } + + case 'h': + if (IC_GKEY(ic) != 1) { + IC_GKEY(ic) = 1 + newgraph = YES + } + + case 'i': + if (IC_GKEY(ic) != 2) { + IC_GKEY(ic) = 2 + newgraph = YES + } + + case 'j': + if (IC_GKEY(ic) != 3) { + IC_GKEY(ic) = 3 + newgraph = YES + } + + case 'k': + if (IC_GKEY(ic) != 4) { + IC_GKEY(ic) = 4 + newgraph = YES + } + + case 'l': + if (IC_GKEY(ic) != 5) { + IC_GKEY(ic) = 5 + newgraph = YES + } + + case 'o': # Set overplot flag + IC_OVERPLOT(ic) = YES + + case 'r': # Redraw the graph + newgraph = YES + + case 'u': # Undelete data points. + call icg_undeleted (ic, gp, gt, cv, Memd[x], Memd[y], + Memd[wts], Memd[userwts], npts, wx, wy) + + case 'w': # Window graph + call gt_window (gt, gp, cursor, newgraph) + + case 'x': # Reset the value of the x point. + i = icg_nearestd (ic, gp, gt, cv, Memd[x], Memd[y], npts, wx, + wy) + + if (i != 0) { + call printf ("Enter new x (%g): ") + call pargd (Memd[x+i-1]) + call flush (STDOUT) + if (scan() != EOF) { + call gargd (x1) + if (nscan() == 1) { + if (!IS_INDEF (x1)) { + oldx = Memd[x+i-1] + oldy = Memd[y+i-1] + Memd[x+i-1] = x1 + call hd_redraw (gp, oldx, oldy, x1, oldy) + IC_NEWX(ic) = YES + } + } + } + } + + case 'y': # Reset the value of the y point. + i = icg_nearestd (ic, gp, gt, cv, Memd[x], Memd[y], npts, wx, + wy) + + if (i != 0) { + call printf ("Enter new y (%g): ") + call pargd (Memd[y+i-1]) + call flush (STDOUT) + if (scan() != EOF) { + call gargd (x1) + if (nscan() == 1) { + if (!IS_INDEF (x1)) { + oldx = Memd[x+i-1] + oldy = Memd[y+i-1] + Memd[y+i-1] = x1 + call hd_redraw (gp, oldx, oldy, oldx, x1) + IC_NEWY(ic) = YES + } + } + } + } + + case 'z': # Reset the weight value of the nearest point + i = icg_nearestd (ic, gp, gt, cv, Memd[x], Memd[y], npts, wx, + wy) + + if (i != 0) { + call printf ("Enter new weight (%g): ") + call pargd (Memd[wts+i-1]) + call flush (STDOUT) + if (scan() != EOF) { + call gargd (x1) + if (nscan() == 1) { + if (!IS_INDEF (x1)) { + Memd[wts+i-1] = x1 + IC_NEWWTS(ic) = YES + } + } + } + } + + default: # Let the user decide on any other keys. + call icg_user (ic, gp, gt, cv, wx, wy, wcs, key, cmd) + } + + # Redraw the graph if necessary. +redraw_ if (newgraph == YES) { + if (IC_AXES(ic, IC_GKEY(ic), 1) != axes[1]) { + axes[1] = IC_AXES(ic, IC_GKEY(ic), 1) + call gt_setr (gt, GTXMIN, INDEF) + call gt_setr (gt, GTXMAX, INDEF) + } + if (IC_AXES(ic, IC_GKEY(ic), 2) != axes[2]) { + axes[2] = IC_AXES(ic, IC_GKEY(ic), 2) + call gt_setr (gt, GTYMIN, INDEF) + call gt_setr (gt, GTYMAX, INDEF) + } + + # Overplot with a different line type + if (IC_OVERPLOT(ic) == YES) + linetype = min ((gstati (gp, G_PLTYPE) + 1), 4) + else + linetype = GL_SOLID + call gseti (gp, G_PLTYPE, linetype) + + call hdic_ebw (ic, Memd[den], Memd[x], Memd[sdev], Memd[ebw], + npts) + + call icg_graphd (ic, gp, gt, cv, Memd[x], Memd[y], Memd[wts], + Memd[ebw], npts) + + newgraph = NO + } + } until (clgcur (cursor, wx, wy, wcs, key, cmd, SZ_LINE) == EOF) + + call mfree (x, TY_DOUBLE) + call mfree (y, TY_DOUBLE) + call mfree (wts, TY_DOUBLE) + call mfree (userwts, TY_DOUBLE) + call mfree (den, TY_DOUBLE) +end diff --git a/noao/imred/dtoi/hdicfit/hdicggraph.x b/noao/imred/dtoi/hdicfit/hdicggraph.x new file mode 100644 index 00000000..dcd9ade3 --- /dev/null +++ b/noao/imred/dtoi/hdicfit/hdicggraph.x @@ -0,0 +1,329 @@ +include <mach.h> +include <gset.h> +include <pkg/gtools.h> +include "hdicfit.h" + +define MSIZE 2. # Mark size +define SZ_TPLOT 7 # String length for plot type +define EB_WTS 10 +define EB_SDEV 11 +define MAX_SZMARKER 4 + +# ICG_GRAPH -- Graph data and fit. + +procedure icg_graphd (ic, gp, gt, cv, x, y, wts, ebw, npts) + +pointer ic # ICFIT pointer +pointer gp # GIO pointer +pointer gt # GTOOLS pointers +pointer cv # Curfit pointer +double x[npts] # Independent variable +double y[npts] # Dependent variable +double wts[npts] # Weights +double ebw[npts] # Half the error bar width +int npts # Number of points + +pointer xout, yout +int nvalues +errchk malloc + +begin + call malloc (xout, npts, TY_DOUBLE) + call malloc (yout, npts, TY_DOUBLE) + + call icg_axesd (ic, gt, cv, 1, x, y, Memd[xout], npts) + call icg_axesd (ic, gt, cv, 2, x, y, Memd[yout], npts) + call icg_paramsd (ic, cv, x, y, wts, npts, gt) + + call icg_g1d (ic, gp, gt, Memd[xout], Memd[yout], wts, ebw, npts) + + if (npts != IC_NFIT(ic)) { + if ((abs (IC_NAVERAGE(ic)) > 1) || (IC_NREJECT(ic) > 0)) { + call realloc (xout, IC_NFIT(ic), TY_DOUBLE) + call realloc (yout, IC_NFIT(ic), TY_DOUBLE) + call icg_axesd (ic, gt, cv, 1, Memd[IC_XFIT(ic)], + Memd[IC_YFIT(ic)], Memd[xout], IC_NFIT(ic)) + call icg_axesd (ic, gt, cv, 2, Memd[IC_XFIT(ic)], + Memd[IC_YFIT(ic)], Memd[yout], IC_NFIT(ic)) + call icg_g2d (ic, gp, gt, Memd[xout], Memd[yout], + IC_NFIT(ic)) + } + + } else if (IC_NREJECT(ic) > 0) + call icg_g2d (ic, gp, gt, Memd[xout], Memd[yout], npts) + + nvalues = NVALS_FIT + call icg_gfd (ic, gp, gt, cv, nvalues) + + # Mark the the sample regions. + call icg_sampled (ic, gp, gt, x, npts, 1) + + call mfree (xout, TY_DOUBLE) + call mfree (yout, TY_DOUBLE) +end + + +# ICG_G1D -- Plot data vector. + +procedure icg_g1d (ic, gp, gt, x, y, wts, ebw, npts) + +pointer ic # Pointer to ic structure +pointer gp # Pointer to graphics stream +pointer gt # Pointer to gtools structure +double x[npts] # Array of independent variables +double y[npts] # Array of dependent variables +double wts[npts] # Array of weights +double ebw[npts] # Error bar half width in WCS (positive density) +int npts # Number of points + +pointer sp, xr, yr, sz, szmk, gt1, gt2 +char tplot[SZ_TPLOT] +int i, xaxis, yaxis, symbols[3], markplot +real size, rmin, rmax, big +bool fp_equald(), streq(), fp_equalr() +int strncmp() +data symbols/GM_PLUS, GM_BOX, GM_CIRCLE/ +include "hdic.com" + +begin + call smark (sp) + call salloc (xr, npts, TY_REAL) + call salloc (yr, npts, TY_REAL) + call salloc (sz, npts, TY_REAL) + call salloc (szmk, npts, TY_REAL) + + call achtdr (x, Memr[xr], npts) + call achtdr (y, Memr[yr], npts) + + xaxis = (IC_AXES (ic, IC_GKEY(ic), 1)) + yaxis = (IC_AXES (ic, IC_GKEY(ic), 2)) + + # Set up gtools structure for deleted (gt1) and added (gt2) points + call gt_copy (gt, gt1) + call gt_setr (gt1, GTXSIZE, MSIZE) + call gt_setr (gt1, GTYSIZE, MSIZE) + call gt_sets (gt1, GTMARK, "cross") + + call gt_copy (gt, gt2) + call gt_setr (gt2, GTXSIZE, MSIZE) + call gt_setr (gt2, GTYSIZE, MSIZE) + call gt_sets (gt2, GTMARK, "circle") + + markplot = NO + call gt_gets (gt, GTTYPE, tplot, SZ_TPLOT) + if (strncmp (tplot, "mark", 4) == 0) + markplot = YES + + if (IC_OVERPLOT(ic) == NO) { + # Start a new plot + call gclear (gp) + + # Set the graph scale and axes + call gascale (gp, Memr[xr], npts, 1) + call gascale (gp, Memr[yr], npts, 2) + + # If plotting HD curve, set wy2 to maxden, which may have + # been updated if a new endpoint was added. + + if (xaxis == 'y' && yaxis == 'u') + call gswind (gp, INDEF, INDEF, INDEF, real (maxden)) + + call gt_swind (gp, gt) + call gt_labax (gp, gt) + } + + # Calculate size of markers if error bars are being used. If the + # weights are being used as the marker size, they are first scaled + # between 1.0 and the maximum marker size. + + call gt_gets (gt, GTMARK, tplot, SZ_TPLOT) + if (streq (tplot, "hebar") || streq (tplot, "vebar")) { + if (IC_EBARS(ic) == EB_WTS) { + call achtdr (wts, Memr[sz], npts) + call alimr (Memr[sz], npts, rmin, rmax) + if (fp_equalr (rmin, rmax)) + call amovr (Memr[sz], Memr[szmk], npts) + else { + big = real (MAX_SZMARKER) + call amapr (Memr[sz], Memr[szmk], npts,rmin,rmax,1.0, big) + } + } else { + call achtdr (ebw, Memr[sz], npts) + call hd_szmk (gp, gt, xaxis, yaxis, tplot, Memr[sz], + Memr[szmk], npts) + } + } else + call amovkr (MSIZE, Memr[szmk], npts) + + do i = 1, npts { + # Check for deleted point + if (fp_equald (wts[i], 0.0D0)) + call gt_plot (gp, gt1, Memr[xr+i-1], Memr[yr+i-1], 1) + + # Check for added point + else if (fp_equald (ebw[i], ADDED_PT)) + call gt_plot (gp, gt2, Memr[xr+i-1], Memr[yr+i-1], 1) + + else { + size = Memr[szmk+i-1] + call gt_setr (gt, GTXSIZE, size) + call gt_setr (gt, GTYSIZE, size) + call gt_plot (gp, gt, Memr[xr+i-1], Memr[yr+i-1], 1) + } + } + + IC_OVERPLOT(ic) = NO + call sfree (sp) +end + + +# HD_SZMK -- Calculate size of error bar markers. This procedure is +# called when the marker type is hebar or vebar and the marker size is +# to be the error in the point, not its weight in the fit. + +procedure hd_szmk (gp, gt, xaxis, yaxis, mark, insz, outsz, npts) + +pointer gp # Pointer to gio structure +pointer gt # Pointer to gtools structure +int xaxis, yaxis # Codes for x and y axis types +char mark[SZ_TPLOT] # Type of marker to use +real insz[npts] # Standard deviations in WCS units +real outsz[npts] # Output size array in NDC units +int npts # Number of points in arrays + +int gt_geti() +char tplot[SZ_TPLOT] +real dx, dy +bool streq() + +begin + # Check validity of axis types + if (xaxis != 'x' && xaxis != 'u' && yaxis != 'x' && yaxis != 'u') { + call eprintf ("Choose graph type with axes 'u' or 'x'; ") + call eprintf ("Using marker size 2.0\n") + call amovkr (MSIZE, outsz, npts) + return + } + + call gt_gets (gt, GTMARK, tplot, SZ_TPLOT) + if (streq (tplot, "hebar")) { + if (yaxis == 'x' || yaxis == 'u') { + call gt_sets (gt, GTMARK, "vebar") + call eprintf ("Marker switched to vebar\n") + # call flush (STDOUT) + } + } else if (streq (tplot, "vebar")) { + if (xaxis == 'x' || xaxis == 'u') { + call gt_sets (gt, GTMARK, "hebar") + call eprintf ("Marker switched to hebar\n") + # call flush (STDOUT) + } + } + + # Need to scale standard deviation from density to NDC units + call ggscale (gp, 0.0, 0.0, dx, dy) + if (gt_geti (gt, GTTRANSPOSE) == NO) + call adivkr (insz, dx, outsz, npts) + else + call adivkr (insz, dy, outsz, npts) +end + + +# ICG_G2D -- Show sample range and rejected data on plot. + +procedure icg_g2d (ic, gp, gt, x, y, npts) + +pointer ic # ICFIT pointer +pointer gp # GIO pointer +pointer gt # GTOOLS pointer +double x[npts], y[npts] # Data points +int npts # Number of data points + +int i +pointer sp, xr, yr, gt1 + +begin + call smark (sp) + call salloc (xr, npts, TY_REAL) + call salloc (yr, npts, TY_REAL) + call achtdr (x, Memr[xr], npts) + call achtdr (y, Memr[yr], npts) + + call gt_copy (gt, gt1) + call gt_sets (gt1, GTTYPE, "mark") + + # Mark the sample points. + + if (abs (IC_NAVERAGE(ic)) > 1) { + call gt_sets (gt1, GTMARK, "plus") + call gt_setr (gt1, GTXSIZE, real (-abs (IC_NAVERAGE(ic)))) + call gt_setr (gt1, GTYSIZE, 1.) + call gt_plot (gp, gt1, Memr[xr], Memr[yr], npts) + } + + # Mark the rejected points. + + if (IC_NREJECT(ic) > 0) { + call gt_sets (gt1, GTMARK, "diamond") + call gt_setr (gt1, GTXSIZE, MSIZE) + call gt_setr (gt1, GTYSIZE, MSIZE) + do i = 1, npts { + if (Memi[IC_REJPTS(ic)+i-1] == YES) + call gt_plot (gp, gt1, Memr[xr+i-1], Memr[yr+i-1], 1) + } + } + + call gt_free (gt1) + call sfree (sp) +end + + +# ICG_GF[RD] -- Overplot the fit on a graph of the data points. A vector +# is constructed and evaluated, then plotted. + +procedure icg_gfd (ic, gp, gt, cv, npts) + +pointer ic # ICFIT pointer +pointer gp # GIO pointer +pointer gt # GTOOL pointer +pointer cv # CURFIT pointer +int npts # Number of points to plot + +pointer sp, xr, yr, x, y, xo, yo, gt1 +int gstati() + +begin + call smark (sp) + + if (IC_FITERROR(ic) == YES) + return + + call salloc (xr, npts, TY_REAL) + call salloc (yr, npts, TY_REAL) + call salloc (x, npts, TY_DOUBLE) + call salloc (y, npts, TY_DOUBLE) + call salloc (xo, npts, TY_DOUBLE) + call salloc (yo, npts, TY_DOUBLE) + + # Calculate big vector of independent variable values. Note value + # of npts can change with this call. + call hdic_gvec (ic, Memd[x], npts, IC_TRANSFORM(ic)) + + # Calculate vector of fit values. + call dcvvector (cv, Memd[x], Memd[y], npts) + + # Convert to user function or transpose axes. Change type to reals + # for plotting. + call icg_axesd (ic, gt, cv, 1, Memd[x], Memd[y], Memd[xo], npts) + call icg_axesd (ic, gt, cv, 2, Memd[x], Memd[y], Memd[yo], npts) + call achtdr (Memd[xo], Memr[xr], npts) + call achtdr (Memd[yo], Memr[yr], npts) + + call gt_copy (gt, gt1) + call gt_sets (gt1, GTTYPE, "line") + call gt_seti (gt1, GTLINE, gstati (gp, G_PLTYPE)) + call gt_plot (gp, gt1, Memr[xr], Memr[yr], npts) + call gt_free (gt1) + + call sfree (sp) +end diff --git a/noao/imred/dtoi/hdicfit/hdicgnearest.x b/noao/imred/dtoi/hdicfit/hdicgnearest.x new file mode 100644 index 00000000..8d556273 --- /dev/null +++ b/noao/imred/dtoi/hdicfit/hdicgnearest.x @@ -0,0 +1,72 @@ +include <mach.h> +include <pkg/gtools.h> + +# ICG_NEAREST -- Find the nearest point to the cursor and return the index. +# The nearest point to the cursor in NDC coordinates is determined. +# The cursor is moved to the nearest point selected. + +int procedure icg_nearestd (ic, gp, gt, cv, x, y, npts, wx, wy) + +pointer ic # ICFIT pointer +pointer gp # GIO pointer +pointer gt # GTOOLS pointer +pointer cv # CURFIT pointer +double x[npts], y[npts] # Data points +int npts # Number of points +real wx, wy # Cursor position + +int pt +pointer sp, xout, yout +int icg_nd(), gt_geti() + +begin + call smark (sp) + call salloc (xout, npts, TY_DOUBLE) + call salloc (yout, npts, TY_DOUBLE) + + call icg_axesd (ic, gt, cv, 1, x, y, Memd[xout], npts) + call icg_axesd (ic, gt, cv, 2, x, y, Memd[yout], npts) + + if (gt_geti (gt, GTTRANSPOSE) == NO) + pt = icg_nd (gp, Memd[xout], Memd[yout], npts, wx, wy) + else + pt = icg_nd (gp, Memd[yout], Memd[xout], npts, wy, wx) + call sfree (sp) + + return (pt) +end + + +# ICG_ND -- ?? + +int procedure icg_nd (gp, x, y, npts, wx, wy) + +pointer gp # GIO pointer +double x[npts], y[npts] # Data points +int npts # Number of points +real wx, wy # Cursor position + +int i, j +real x0, y0, r2, r2min + +begin + # Transform world cursor coordinates to NDC. + call gctran (gp, wx, wy, wx, wy, 1, 0) + + # Search for nearest point. + r2min = MAX_REAL + do i = 1, npts { + call gctran (gp, real (x[i]), real (y[i]), x0, y0, 1, 0) + r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2 + if (r2 < r2min) { + r2min = r2 + j = i + } + } + + # Move the cursor to the selected point and return the index. + if (j != 0) + call gscur (gp, real (x[j]), real (y[j])) + + return (j) +end diff --git a/noao/imred/dtoi/hdicfit/hdicgparams.x b/noao/imred/dtoi/hdicfit/hdicgparams.x new file mode 100644 index 00000000..e502fba1 --- /dev/null +++ b/noao/imred/dtoi/hdicfit/hdicgparams.x @@ -0,0 +1,94 @@ +include <pkg/gtools.h> +include "hdicfit.h" + +# ICG_PARAMS -- Set parameter string. + +procedure icg_paramsd (ic, cv, x, y, wts, npts, gt) + +pointer ic # ICFIT pointer +pointer cv # Curfit pointer +double x[ARB] # Ordinates +double y[ARB] # Abscissas +double wts[ARB] # Weights +int npts # Number of data points +pointer gt # GTOOLS pointer + +double rms +int i, n, deleted +pointer sp, fit, wts1, str, params +double ic_rmsd() + +begin + call smark (sp) + + if (npts == IC_NFIT(ic)) { + # Allocate memory for the fit. + n = npts + call salloc (fit, n, TY_DOUBLE) + call salloc (wts1, n, TY_DOUBLE) + + # Eliminate rejected points and count deleted points. + call amovd (wts, Memd[wts1], n) + if (IC_NREJECT(ic) > 0) { + do i = 1, npts { + if (Memi[IC_REJPTS(ic)+i-1] == YES) + Memd[wts1+i-1] = 0. + } + } + deleted = 0 + do i = 1, n { + if (wts[i] == 0.) + deleted = deleted + 1 + } + + # Set the fit and compute the RMS error. + call dcvvector (cv, x, Memd[fit], n) + rms = ic_rmsd (x, y, Memd[fit], Memd[wts1], n) + + } else { + # Allocate memory for the fit. + n = IC_NFIT(ic) + call salloc (fit, n, TY_DOUBLE) + call salloc (wts1, n, TY_DOUBLE) + + # Eliminate rejected points and count deleted points. + call amovd (Memd[IC_WTSFIT(ic)], Memd[wts1], n) + if (IC_NREJECT(ic) > 0) { + do i = 1, npts { + if (Memi[IC_REJPTS(ic)+i-1] == YES) + Memd[wts1+i-1] = 0. + } + } + deleted = 0 + do i = 1, n { + if (wts[i] == 0.) + deleted = deleted + 1 + } + + # Set the fit and compute the rms error. + call dcvvector (cv, Memd[IC_XFIT(ic)], Memd[fit], n) + rms = ic_rmsd (Memd[IC_XFIT(ic)], Memd[IC_YFIT(ic)], Memd[fit], + Memd[wts1], n) + } + + # Print the parameters and errors. + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (params, 2*SZ_LINE, TY_CHAR) + + call sprintf (Memc[str],SZ_LINE, "function=%s, order=%d, transform=%s") + call ic_gstr (ic, "function", Memc[params], 2*SZ_LINE) + call pargstr (Memc[params]) + call pargi (IC_ORDER(ic)) + call ic_gstr (ic, "transform", Memc[params], SZ_LINE) + call pargstr (Memc[params]) + call sprintf (Memc[params], 2*SZ_LINE, + "%s\nfog=%.5f, total=%d, deleted=%d, RMS=%7.4g") + call pargstr (Memc[str]) + call pargr (IC_FOG(ic)) + call pargi (npts) + call pargi (deleted) + call pargd (rms) + call gt_sets (gt, GTPARAMS, Memc[params]) + + call sfree (sp) +end diff --git a/noao/imred/dtoi/hdicfit/hdicgredraw.x b/noao/imred/dtoi/hdicfit/hdicgredraw.x new file mode 100644 index 00000000..d42a0740 --- /dev/null +++ b/noao/imred/dtoi/hdicfit/hdicgredraw.x @@ -0,0 +1,22 @@ +include <gset.h> + +define MSIZE 2.0 + +# HD_REDRAW -- Redraw a data point. The old position marker is erased, +# the new marker drawn, and the cursor moved to the new location. + +procedure hd_redraw (gp, oldx, oldy, newx, newy) + +pointer gp # Pointer to graphics stream +real oldx # Old x coordinate +real oldy # Old y coordinate +real newx # New x coordinate +real newy # New y coordinate + +begin + call gseti (gp, G_PMLTYPE, GL_CLEAR) + call gmark (gp, oldx, oldy, GM_PLUS, MSIZE, MSIZE) + call gseti (gp, G_PMLTYPE, GL_SOLID) + call gmark (gp, newx, newy, GM_PLUS, MSIZE, MSIZE) + call gscur (gp, newx, newy) +end diff --git a/noao/imred/dtoi/hdicfit/hdicgsample.x b/noao/imred/dtoi/hdicfit/hdicgsample.x new file mode 100644 index 00000000..12eef158 --- /dev/null +++ b/noao/imred/dtoi/hdicfit/hdicgsample.x @@ -0,0 +1,84 @@ +include <gset.h> +include <pkg/rg.h> +include <pkg/gtools.h> +include "hdicfit.h" + +# ICG_SAMPLE -- Mark sample. + +procedure icg_sampled (ic, gp, gt, x, npts, pltype) + +pointer ic # ICFIT pointer +pointer gp # GIO pointer +pointer gt # GTOOLS pointer +double x[npts] # Ordinates of graph +int npts # Number of data points +int pltype # Plot line type + +pointer rg +int i, axis, pltype1 +real xl, xr, yb, yt, dy +real x1, x2, y1, y2, y3 + +int gstati(), stridxs(), gt_geti() +pointer rg_xrangesd() + +begin + if (stridxs ("*", Memc[IC_SAMPLE(ic)]) > 0) + return + + # Find axis along which the independent data is plotted. + if (IC_AXES(ic,IC_GKEY(ic),1) == 'x') + axis = 1 + else if (IC_AXES(ic,IC_GKEY(ic),2) == 'x') + axis = 2 + else + return + + if (gt_geti (gt, GTTRANSPOSE) == YES) + axis = mod (axis, 2) + 1 + + pltype1 = gstati (gp, G_PLTYPE) + call gseti (gp, G_PLTYPE, pltype) + rg = rg_xrangesd (Memc[IC_SAMPLE(ic)], x, npts) + + switch (axis) { + case 1: + call ggwind (gp, xl, xr, yb, yt) + + dy = yt - yb + y1 = yb + dy / 100 + y2 = y1 + dy / 20 + y3 = (y1 + y2) / 2 + + do i = 1, RG_NRGS(rg) { + x1 = x[RG_X1(rg, i)] + x2 = x[RG_X2(rg, i)] + if ((x1 > xl) && (x1 < xr)) + call gline (gp, x1, y1, x1, y2) + if ((x2 > xl) && (x2 < xr)) + call gline (gp, x2, y1, x2, y2) + call gline (gp, x1, y3, x2, y3) + } + case 2: + call ggwind (gp, yb, yt, xl, xr) + + dy = yt - yb + y1 = yb + dy / 100 + y2 = y1 + dy / 20 + y3 = (y1 + y2) / 2 + + do i = 1, RG_NRGS(rg) { + x1 = x[RG_X1(rg, i)] + x2 = x[RG_X2(rg, i)] + if ((x1 > xl) && (x1 < xr)) + call gline (gp, y1, x1, y2, x1) + if ((x2 > xl) && (x2 < xr)) + call gline (gp, y1, x2, y2, x2) + call gline (gp, y3, x1, y3, x2) + } + } + + call gseti (gp, G_PLTYPE, pltype1) + call rg_free (rg) +end + diff --git a/noao/imred/dtoi/hdicfit/hdicguaxes.x b/noao/imred/dtoi/hdicfit/hdicguaxes.x new file mode 100644 index 00000000..f8199c87 --- /dev/null +++ b/noao/imred/dtoi/hdicfit/hdicguaxes.x @@ -0,0 +1,38 @@ +include "hdicfit.h" + +# ICG_UAXIS -- Set user axis. + +procedure icg_uaxesd (ic, key, cv, x, y, z, npts, label, units, maxchars) + +pointer ic # Pointer to ic structure +int key # Key for axes +pointer cv # CURFIT pointer +double x[npts] # Independent variable +double y[npts] # Dependent variable +double z[npts] # Output values +int npts # Number of points +char label[maxchars] # Axis label +char units[maxchars] # Units for axis +int maxchars # Maximum chars in label + +int offset +double fog +real ic_getr() +include "hdic.com" + +begin + # Axis type 'u' returns the untransformed independent variable + # in the z array. That is, the original density values after + # subtracting the current fog value. Some density values could be + # below fog were excluded from the transformed vector. + + call strcpy ("Density above fog", label, maxchars) + fog = double (ic_getr (ic, "fog")) + + if (npts == nraw) + call asubkd (Memd[den], fog, z, npts) + else { + offset = big_den + (NVALS_FIT - npts) + call asubkd (Memd[offset], fog, z, npts) + } +end diff --git a/noao/imred/dtoi/hdicfit/hdicgundel.x b/noao/imred/dtoi/hdicfit/hdicgundel.x new file mode 100644 index 00000000..9a7c68a4 --- /dev/null +++ b/noao/imred/dtoi/hdicfit/hdicgundel.x @@ -0,0 +1,87 @@ +include <gset.h> +include <mach.h> +include <pkg/gtools.h> +include "hdicfit.h" + +define MSIZE 2. # Mark size + +# ICG_UNDELETE -- Undelete data point nearest the cursor. +# The nearest point to the cursor in NDC coordinates is determined. + +procedure icg_undeleted (ic, gp, gt, cv, x, y, wts, userwts, npts, wx, wy) + +pointer ic # ICFIT pointer +pointer gp # GIO pointer +pointer gt # GTOOLS pointer +pointer cv # CURFIT pointer +double x[npts], y[npts] # Data points +double wts[npts], userwts[npts] # Weight arrays +int npts # Number of points +real wx, wy # Position to be nearest + +pointer sp, xout, yout +int gt_geti() + +begin + call smark (sp) + call salloc (xout, npts, TY_DOUBLE) + call salloc (yout, npts, TY_DOUBLE) + + call icg_axesd (ic, gt, cv, 1, x, y, Memd[xout], npts) + call icg_axesd (ic, gt, cv, 2, x, y, Memd[yout], npts) + if (gt_geti (gt, GTTRANSPOSE) == NO) { + call icg_u1d (ic, gp, Memd[xout], Memd[yout], wts, userwts, + npts, wx, wy) + } else { + call icg_u1d (ic, gp, Memd[yout], Memd[xout], wts, userwts, + npts, wy, wx) + } + + call sfree (sp) +end + + +# ICG_U1D -- Do the actual undelete. + +procedure icg_u1d (ic, gp, x, y, wts, userwts, npts, wx, wy) + +pointer ic # ICFIT pointer +pointer gp # GIO pointer +double x[npts], y[npts] # Data points +double wts[npts], userwts[npts] # Weight arrays +int npts # Number of points +real wx, wy # Position to be nearest + +int i, j +real x0, y0, r2, r2min + +begin + # Transform world cursor coordinates to NDC. + call gctran (gp, wx, wy, wx, wy, 1, 0) + + # Search for nearest point to a point with zero weight. + r2min = MAX_REAL + do i = 1, npts { + if (wts[i] != 0.) + next + + call gctran (gp, real (x[i]), real (y[i]), x0, y0, 1, 0) + + r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2 + if (r2 < r2min) { + r2min = r2 + j = i + } + } + + # Unmark the deleted point and reset the weight. + if (j != 0) { + call gscur (gp, real (x[j]), real (y[j])) + call gseti (gp, G_PMLTYPE, GL_CLEAR) + call gmark (gp, real (x[j]), real (y[j]), GM_CROSS, MSIZE, MSIZE) + call gseti (gp, G_PMLTYPE, GL_SOLID) + call gmark (gp, real (x[j]), real (y[j]), GM_PLUS, MSIZE, MSIZE) + wts[j] = userwts[j] + IC_NEWWTS(ic) = YES + } +end diff --git a/noao/imred/dtoi/hdicfit/hdicguser.x b/noao/imred/dtoi/hdicfit/hdicguser.x new file mode 100644 index 00000000..5e070537 --- /dev/null +++ b/noao/imred/dtoi/hdicfit/hdicguser.x @@ -0,0 +1,17 @@ +# ICG_USER -- User default action + +procedure icg_user (ic, gp, gt, cv, wx, wy, wcs, key, cmd) + +pointer ic # ICFIT pointer +pointer gp # GIO pointer +pointer gt # GTOOLS pointer +pointer cv # CURFIT pointer +real wx, wy # Cursor positions +int wcs # GIO WCS +int key # Cursor key +char cmd[ARB] # Cursor command + +begin + # Ring bell + call printf ("\07") +end diff --git a/noao/imred/dtoi/hdicfit/hdicgvec.x b/noao/imred/dtoi/hdicfit/hdicgvec.x new file mode 100644 index 00000000..81909cb8 --- /dev/null +++ b/noao/imred/dtoi/hdicfit/hdicgvec.x @@ -0,0 +1,74 @@ +include <mach.h> +include "hdicfit.h" + +# HDIC_GVEC -- Get a vector of data to be plotted. The raw density +# values are stored in common. From these values, all possible transforms +# can be generated. + +procedure hdic_gvec (ic, xout, npts, transform) + +pointer ic # Pointer to ic structure +double xout[npts] # Output vector +int npts # Desired npoints to output - possibly changed on output +int transform # Integer code for desired type of transform + +pointer sp, bigdenaf +int i, j, npts_desired +double fog, dval +real ic_getr() +include "hdic.com" + +begin + npts_desired = npts + if (npts_desired != NVALS_FIT) + call error (0, "hdicgvec: nvals != NVALS_FIT") + + call smark (sp) + call salloc (bigdenaf, npts, TY_DOUBLE) + + j = 1 + + fog = double (ic_getr (ic, "fog")) + call asubkd (Memd[big_den], fog, Memd[bigdenaf], npts) + + switch (transform) { + case HD_NONE: + call amovd (Memd[bigdenaf], xout, NVALS_FIT) + + case HD_LOGO: + do i = 1, npts_desired { + dval = Memd[bigdenaf+i-1] + if (dval < 0.0D0) + npts = npts - 1 + else { + xout[j] = log10 ((10.**dval) - 1.0) + j = j + 1 + } + } + + case HD_K75: + do i = 1, npts_desired { + dval = Memd[bigdenaf+i-1] + if (dval < 0.0D0) + npts = npts - 1 + else { + xout[j] = dval + 0.75 * log10 (1.0 - (10.** (-dval))) + j = j + 1 + } + } + + case HD_K50: + do i = 1, npts_desired { + dval = Memd[bigdenaf+i-1] + if (dval < 0.0D0) + npts = npts - 1 + else { + xout[j] = dval + 0.50 * log10 (1.0 - (10.** (-dval))) + j = j + 1 + } + } + + } + + call sfree (sp) +end diff --git a/noao/imred/dtoi/hdicfit/hdicinit.x b/noao/imred/dtoi/hdicfit/hdicinit.x new file mode 100644 index 00000000..d9730051 --- /dev/null +++ b/noao/imred/dtoi/hdicfit/hdicinit.x @@ -0,0 +1,60 @@ +include <mach.h> +include "hdicfit.h" + +# HDIC_INIT -- Initialize hdfit/icgfit interface. + +procedure hdic_init (density, nvalues, dmax) + +double density[nvalues] # Reference density above fog values +int nvalues # Number of values in sample +double dmax # Maximum possible density + +int i +pointer sp, base +double xxmax, xxmin, delta_den +bool fp_equald() +include "hdic.com" +errchk malloc + +begin + call smark (sp) + + if (den == NULL || big_den == NULL) { + call malloc (den, nvalues, TY_DOUBLE) + call malloc (big_den, NVALS_FIT, TY_DOUBLE) + } else if (nvalues != nraw) + call realloc (den, nvalues, TY_DOUBLE) + + nraw = nvalues + call salloc (base, NVALS_FIT, TY_DOUBLE) + + # Copy density array to pointer location + call amovd (density, Memd[den], nraw) + + # Calculate big vector of density values. The points are spaced + # linear in log space, to yield adequate spacing at low density values. + + call alimd (density, nraw, xxmin, xxmax) + + # Put user value for maximum density in common block if it is valid. + if (! fp_equald (dmax, 0.0D0)) + maxden = dmax + + # Make big_den go all the way to maxden, not just xxmax. Make sure + # the value of xxmin won't cause the log function to blow up. + + if (xxmin > 0.0D0) + ; + else + xxmin = 2.0 * EPSILOND + + delta_den = (log10 (maxden) - log10 (xxmin)) / double (NVALS_FIT - 1) + + do i = 1, NVALS_FIT + Memd[big_den+i-1] = log10 (xxmin) + double (i-1) * delta_den + + call amovkd (10.0D0, Memd[base], NVALS_FIT) + call aexpd (Memd[base], Memd[big_den], Memd[big_den], NVALS_FIT) + + call sfree (sp) +end diff --git a/noao/imred/dtoi/hdicfit/hdicparams.x b/noao/imred/dtoi/hdicfit/hdicparams.x new file mode 100644 index 00000000..1840dccd --- /dev/null +++ b/noao/imred/dtoi/hdicfit/hdicparams.x @@ -0,0 +1,323 @@ +include "hdicfit.h" + +define FUNCTIONS "|chebyshev|legendre|spline3|spline1|power|" +define TRANSFORMS "|none|logopacitance|k50|k75|" + +# IC_OPEN -- Open ICFIT parameter structure. + +procedure ic_open (ic) + +pointer ic # ICFIT pointer +errchk malloc + +begin + # Allocate memory for the package parameter structure. + call malloc (ic, IC_LENSTRUCT, TY_STRUCT) + call malloc (IC_SAMPLE(ic), SZ_LINE, TY_CHAR) + call malloc (IC_LABELS(ic,1), SZ_LINE, TY_CHAR) + call malloc (IC_LABELS(ic,2), SZ_LINE, TY_CHAR) + call malloc (IC_UNITS(ic,1), SZ_LINE, TY_CHAR) + call malloc (IC_UNITS(ic,2), SZ_LINE, TY_CHAR) + + # Set defaults. + call ic_pstr (ic, "function", "spline3") + call ic_puti (ic, "order", 1) + call ic_pstr (ic, "sample", "*") + call ic_puti (ic, "naverage", 1) + call ic_puti (ic, "niterate", 0) + call ic_putr (ic, "low", 0.) + call ic_putr (ic, "high", 0.) + call ic_putr (ic, "grow", 0.) + call ic_pstr (ic, "xlabel", "X") + call ic_pstr (ic, "ylabel", "Y") + call ic_pstr (ic, "xunits", "") + call ic_pstr (ic, "yunits", "") + call ic_puti (ic, "key", 1) + call ic_pkey (ic, 1, 'x', 'y') + call ic_pkey (ic, 2, 'y', 'x') + call ic_pkey (ic, 3, 'x', 'r') + call ic_pkey (ic, 4, 'x', 'd') + call ic_pkey (ic, 5, 'x', 'n') + + # Initialize other parameters + IC_RG(ic) = NULL + IC_XFIT(ic) = NULL + IC_YFIT(ic) = NULL + IC_WTSFIT(ic) = NULL + IC_REJPTS(ic) = NULL +end + + +# IC_CLOSER -- Close ICFIT parameter structure. + +procedure ic_closer (ic) + +pointer ic # ICFIT pointer + +begin + # Free memory for the package parameter structure. + call rg_free (IC_RG(ic)) + call mfree (IC_XFIT(ic), TY_REAL) + call mfree (IC_YFIT(ic), TY_REAL) + call mfree (IC_WTSFIT(ic), TY_REAL) + call mfree (IC_REJPTS(ic), TY_INT) + call mfree (IC_SAMPLE(ic), TY_CHAR) + call mfree (IC_LABELS(ic,1), TY_CHAR) + call mfree (IC_LABELS(ic,2), TY_CHAR) + call mfree (IC_UNITS(ic,1), TY_CHAR) + call mfree (IC_UNITS(ic,2), TY_CHAR) + call mfree (ic, TY_STRUCT) +end + + +# IC_CLOSED -- Close ICFIT parameter structure. + +procedure ic_closed (ic) + +pointer ic # ICFIT pointer + +begin + # Free memory for the package parameter structure. + call rg_free (IC_RG(ic)) + call mfree (IC_XFIT(ic), TY_DOUBLE) + call mfree (IC_YFIT(ic), TY_DOUBLE) + call mfree (IC_WTSFIT(ic), TY_DOUBLE) + call mfree (IC_REJPTS(ic), TY_INT) + call mfree (IC_SAMPLE(ic), TY_CHAR) + call mfree (IC_LABELS(ic,1), TY_CHAR) + call mfree (IC_LABELS(ic,2), TY_CHAR) + call mfree (IC_UNITS(ic,1), TY_CHAR) + call mfree (IC_UNITS(ic,2), TY_CHAR) + call mfree (ic, TY_STRUCT) +end + + +# IC_PSTR -- Put string valued parameters. + +procedure ic_pstr (ic, param, str) + +pointer ic # ICFIT pointer +char param[ARB] # Parameter to be put +char str[ARB] # String value + +int i +pointer ptr +int strdic() +bool streq() + +begin + if (streq (param, "sample")) + call strcpy (str, Memc[IC_SAMPLE(ic)], SZ_LINE) + else if (streq (param, "function")) { + call malloc (ptr, SZ_LINE, TY_CHAR) + i = strdic (str, Memc[ptr], SZ_LINE, FUNCTIONS) + if (i > 0) + IC_FUNCTION(ic) = i + call mfree (ptr, TY_CHAR) + } else if (streq (param, "transform")) { + call malloc (ptr, SZ_LINE, TY_CHAR) + i = strdic (str, Memc[ptr], SZ_LINE, TRANSFORMS) + if (i > 0) + IC_TRANSFORM(ic) = i + call mfree (ptr, TY_CHAR) + } else if (streq (param, "xlabel")) + call strcpy (str, Memc[IC_LABELS(ic,1)], SZ_LINE) + else if (streq (param, "ylabel")) + call strcpy (str, Memc[IC_LABELS(ic,2)], SZ_LINE) + else if (streq (param, "xunits")) + call strcpy (str, Memc[IC_UNITS(ic,1)], SZ_LINE) + else if (streq (param, "yunits")) + call strcpy (str, Memc[IC_UNITS(ic,2)], SZ_LINE) + else + call error (0, "ICFIT: Unknown parameter") +end + + +# IC_PUTI -- Put integer valued parameters. + +procedure ic_puti (ic, param, ival) + +pointer ic # ICFIT pointer +char param[ARB] # Parameter to be put +int ival # Integer value + +bool streq() + +begin + if (streq (param, "naverage")) + IC_NAVERAGE(ic) = ival + else if (streq (param, "order")) + IC_ORDER(ic) = ival + else if (streq (param, "niterate")) + IC_NITERATE(ic) = ival + else if (streq (param, "key")) + IC_GKEY(ic) = ival + else if (streq (param, "transform")) + IC_TRANSFORM(ic) = ival + else + call error (0, "ICFIT: Unknown parameter") +end + + +# IC_PKEY -- Put key parameters. +# Note the key types must be integers not characters. + +procedure ic_pkey (ic, key, xaxis, yaxis) + +pointer ic # ICFIT pointer +int key # Key to be defined +int xaxis # X axis type +int yaxis # Y axis type + +begin + IC_AXES(ic, key, 1) = xaxis + IC_AXES(ic, key, 2) = yaxis +end + + +# IC_GKEY -- Get key parameters. + +procedure ic_gkey (ic, key, xaxis, yaxis) + +pointer ic # ICFIT pointer +int key # Key to be gotten +int xaxis # X axis type +int yaxis # Y axis type + +begin + xaxis = IC_AXES(ic, key, 1) + yaxis = IC_AXES(ic, key, 2) +end + + +# IC_PUTR -- Put real valued parameters. + +procedure ic_putr (ic, param, rval) + +pointer ic # ICFIT pointer +char param[ARB] # Parameter to be put +real rval # Real value + +bool streq() + +begin + if (streq (param, "xmin")) + IC_XMIN(ic) = rval + else if (streq (param, "xmax")) + IC_XMAX(ic) = rval + else if (streq (param, "low")) + IC_LOW(ic) = rval + else if (streq (param, "high")) + IC_HIGH(ic) = rval + else if (streq (param, "grow")) + IC_GROW(ic) = rval + else if (streq (param, "fog")) + IC_FOG(ic) = rval + else if (streq (param, "rfog")) + IC_RFOG(ic) = rval + else + call error (0, "ICFIT: Unknown parameter") +end + + +# IC_GSTR -- Get string valued parameters. + +procedure ic_gstr (ic, param, str, maxchars) + +pointer ic # ICFIT pointer +char param[ARB] # Parameter to be put +char str[maxchars] # String value +int maxchars # Maximum number of characters + +bool streq() + +begin + if (streq (param, "sample")) + call strcpy (Memc[IC_SAMPLE(ic)], str, maxchars) + else if (streq (param, "xlabel")) + call strcpy (Memc[IC_LABELS(ic,1)], str, maxchars) + else if (streq (param, "ylabel")) + call strcpy (Memc[IC_LABELS(ic,2)], str, maxchars) + else if (streq (param, "xunits")) + call strcpy (Memc[IC_UNITS(ic,1)], str, maxchars) + else if (streq (param, "yunits")) + call strcpy (Memc[IC_UNITS(ic,2)], str, maxchars) + else if (streq (param, "function")) { + switch (IC_FUNCTION(ic)) { + case 1: + call strcpy ("chebyshev", str, maxchars) + case 2: + call strcpy ("legendre", str, maxchars) + case 3: + call strcpy ("spline3", str, maxchars) + case 4: + call strcpy ("spline1", str, maxchars) + case 5: + call strcpy ("power", str, maxchars) + } + } else if (streq (param, "transform")) { + switch (IC_TRANSFORM(ic)) { + case 1: + call strcpy ("none", str, maxchars) + case 2: + call strcpy ("logopacitance", str, maxchars) + case 3: + call strcpy ("k50", str, maxchars) + case 4: + call strcpy ("k75", str, maxchars) + } + } else + call error (0, "ICFIT: Unknown parameter") +end + + +# IC_GETI -- Get integer valued parameters. + +int procedure ic_geti (ic, param) + +pointer ic # ICFIT pointer +char param[ARB] # Parameter to be gotten + +bool streq() + +begin + if (streq (param, "naverage")) + return (IC_NAVERAGE(ic)) + else if (streq (param, "order")) + return (IC_ORDER(ic)) + else if (streq (param, "niterate")) + return (IC_NITERATE(ic)) + else if (streq (param, "key")) + return (IC_GKEY(ic)) + else if (streq (param, "transform")) + return (IC_TRANSFORM(ic)) + + call error (0, "ICFIT: Unknown parameter") +end + + +# IC_GETR -- Get real valued parameters. + +real procedure ic_getr (ic, param) + +pointer ic # ICFIT pointer +char param[ARB] # Parameter to be put + +bool streq() + +begin + if (streq (param, "xmin")) + return (IC_XMIN(ic)) + else if (streq (param, "xmax")) + return (IC_XMAX(ic)) + else if (streq (param, "low")) + return (IC_LOW(ic)) + else if (streq (param, "high")) + return (IC_HIGH(ic)) + else if (streq (param, "grow")) + return (IC_GROW(ic)) + else if (streq (param, "fog")) + return (IC_FOG(ic)) + + call error (0, "ICFIT: Unknown parameter") +end diff --git a/noao/imred/dtoi/hdicfit/hdicreject.x b/noao/imred/dtoi/hdicfit/hdicreject.x new file mode 100644 index 00000000..82dd093c --- /dev/null +++ b/noao/imred/dtoi/hdicfit/hdicreject.x @@ -0,0 +1,39 @@ +# IC_REJECT -- Reject points with large residuals from the fit. +# +# The sigma of the fit residuals is calculated. The rejection thresholds +# are set at low_reject*sigma and high_reject*sigma. Points outside the +# rejection threshold are rejected from the fit and flagged in the rejpts +# array. Finally, the remaining points are refit. + +procedure ic_rejectd (cv, x, y, w, rejpts, npts, low_reject, high_reject, + niterate, grow, nreject) + +pointer cv # Curve descriptor +double x[npts] # Input ordinates +double y[npts] # Input data values +double w[npts] # Weights +int rejpts[npts] # Points rejected +int npts # Number of input points +real low_reject, high_reject # Rejection threshold +int niterate # Number of rejection iterations +real grow # Rejection radius +int nreject # Number of points rejected + +int i, newreject + +begin + # Initialize rejection. + nreject = 0 + call amovki (NO, rejpts, npts) + + if (niterate <= 0) + return + + # Find deviant points. + do i = 1, niterate { + call ic_deviantd (cv, x, y, w, rejpts, npts, low_reject, + high_reject, grow, YES, nreject, newreject) + if (newreject == 0) + break + } +end diff --git a/noao/imred/dtoi/hdicfit/hdicshow.x b/noao/imred/dtoi/hdicfit/hdicshow.x new file mode 100644 index 00000000..521f04d3 --- /dev/null +++ b/noao/imred/dtoi/hdicfit/hdicshow.x @@ -0,0 +1,52 @@ +include <pkg/gtools.h> +include "hdicfit.h" + +# IC_SHOW -- Show the values of the parameters. + +procedure ic_show (ic, file, gt) + +pointer ic # ICFIT pointer +char file[ARB] # Output file +pointer gt # GTOOLS pointer + +int fd +pointer str +int open() +long clktime() +errchk open, malloc + +begin + fd = open (file, APPEND, TEXT_FILE) + call malloc (str, SZ_LINE, TY_CHAR) + + call cnvtime (clktime(0), Memc[str], SZ_LINE) + call fprintf (fd, "\n# %s\n") + call pargstr (Memc[str]) + + call gt_gets (gt, GTTITLE, Memc[str], SZ_LINE) + call fprintf (fd, "# %s\n") + call pargstr (Memc[str]) + + call gt_gets (gt, GTYUNITS, Memc[str], SZ_LINE) + if (Memc[str] != EOS) { + call fprintf (fd, "fit units = %s\n") + call pargstr (Memc[str]) + } + + call ic_gstr (ic, "function", Memc[str], SZ_LINE) + call fprintf (fd, "function = %s\n") + call pargstr (Memc[str]) + + call fprintf (fd, "order = %d\n") + call pargi (IC_ORDER(ic)) + + call ic_gstr (ic, "transform", Memc[str], SZ_LINE) + call fprintf (fd, "transform = %s\n") + call pargstr (Memc[str]) + + call fprintf (fd, "fog = %g\n") + call pargr (IC_FOG(ic)) + + call mfree (str, TY_CHAR) + call close (fd) +end diff --git a/noao/imred/dtoi/hdicfit/hdicsort.x b/noao/imred/dtoi/hdicfit/hdicsort.x new file mode 100644 index 00000000..16d6e0a6 --- /dev/null +++ b/noao/imred/dtoi/hdicfit/hdicsort.x @@ -0,0 +1,38 @@ +# HDIC_SORT -- sort the log exposure, density and weight information in order +# of increasing density value. The sorting is done is place. The four +# data values are assummed matched on input, that is, exposure[i] matches +# density[i] with weight[i] (and userwts[i]) for all array entries. + +procedure hdic_sort (density, exposure, weights, userwts, whydel, sdev, nvals) + +double density[nvals] # Density array +double exposure[nvals] # Log exposure array +double weights[nvals] # Weights array +double userwts[nvals] # Reference weights array +int whydel[nvals] # Flag array of reasons for deletion +double sdev[nvals] # Array of standard deviations +int nvals # Number of values to sort + +int i, j +double temp +define swap {temp=$1;$1=$2;$2=temp} +int itemp +define iswap {itemp=$1;$1=$2;$2=itemp} + +begin + # Bubble sort - inefficient, but sorting is done infrequently + # an expected small sample size (16 pts typically). + + for (i = nvals; i > 1; i = i - 1) + for (j = 1; j < i; j = j + 1) + if (density [j] > density[j+1]) { + + # Out of order; exchange values + swap (exposure[j], exposure[j+1]) + swap ( density[j], density[j+1]) + swap ( weights[j], weights[j+1]) + swap ( userwts[j], userwts[j+1]) + iswap ( whydel[j], whydel[j+1]) + swap ( sdev[j], sdev[j+1]) + } +end diff --git a/noao/imred/dtoi/hdicfit/hdictrans.x b/noao/imred/dtoi/hdicfit/hdictrans.x new file mode 100644 index 00000000..0ee89037 --- /dev/null +++ b/noao/imred/dtoi/hdicfit/hdictrans.x @@ -0,0 +1,155 @@ +include <mach.h> +include "hdicfit.h" + +# HDIC_TRANSFORM -- Transform density to independent variable of fit. The +# desired transform is stored in the ic structure. A vector of x values +# is returned, as is a possibly modified weights array. The minimum and +# maximum limits of the fit are updated in the ic structure; the labels +# are set also when IC_NEWTRANSFORM = YES. The fog value is subtracted +# from the input density array and the transform performed. + +procedure hdic_transform (ic, density, userwts, xout, wts, whydel, npts) + +pointer ic # Pointer to ic structure +double density[npts] # Array of original density values +double userwts[npts] # Array of original weights values +double xout[npts] # Transformed density above fog (returned) +double wts[npts] # Input weights array +int whydel[npts] # Reason for deletion array +int npts # The number of density points - maybe changed on output + +int i +pointer denaf, sp +double fog, xxmin, xxmax, dval +bool fp_equald() +real ic_getr() +include "hdic.com" + +begin + # Allocate space for density above fog array + call smark (sp) + call salloc (denaf, npts, TY_DOUBLE) + + fog = double (ic_getr (ic, "fog")) + call asubkd (density, fog, Memd[denaf], npts) + + switch (IC_TRANSFORM(ic)) { + case HD_NONE: + do i = 1, npts { + xout[i] = Memd[denaf+i-1] + # In every case, if the point was deleted by the program, + # restore it. + if (whydel[i] == PDELETE) { + wts[i] = userwts[i] + whydel[i] = NDELETE + } + } + + call ic_pstr (ic, "xlabel", "Density above Fog") + xxmin = MIN_DEN - fog + xxmax = maxden + call ic_putr (ic, "xmin", real (xxmin)) + call ic_putr (ic, "xmax", real (xxmax)) + + case HD_LOGO: + call ic_pstr (ic, "xlabel", "Log Opacitance: Log (10**Den - 1)") + xxmin = log10 ((10. ** (MIN_DEN)) - 1.0) + xxmax = log10 ((10. ** (maxden)) - 1.0) + call ic_putr (ic, "xmin", real (xxmin)) + call ic_putr (ic, "xmax", real (xxmax)) + + do i = 1, npts { + dval = Memd[denaf+i-1] + if (dval < 0.0D0 || (fp_equald (dval, 0.0D0))) { + xout[i] = dval + wts[i] = 0.0D0 + whydel[i] = PDELETE + + } else { + xout[i] = log10 ((10. ** (dval)) - 1.0) + + # If point had been deleted, find out why. It affects the + # weights value returned. Only if the point was previously + # deleted by the program, restore it; otherwise, leave it + # alone. + + if (fp_equald (wts[i], 0.0D0)) { + if (whydel[i] == PDELETE) { + wts[i] = userwts[i] + whydel[i] = NDELETE + } + } else + wts[i] = userwts[i] + } + } + + case HD_K75: + call ic_pstr (ic, "xlabel", "Den + 0.75 * Log (1 - (10 ** -Den))") + xxmin = MIN_DEN + 0.75 * log10 (1.0 - (10. ** (-MIN_DEN))) + xxmax = maxden + 0.75 * log10 (1.0 - (10. ** (-maxden))) + call ic_putr (ic, "xmin", real (xxmin)) + call ic_putr (ic, "xmax", real (xxmax)) + + do i = 1, npts { + dval = Memd[denaf+i-1] + if (dval < 0.0D0 || (fp_equald (dval, 0.0D0))) { + xout[i] = dval + wts[i] = 0.0D0 + whydel[i] = PDELETE + + } else { + xout[i] = dval + 0.75 * log10 (1.0 - (10.** (-dval))) + + # If point had been deleted, find out why. It affects the + # weights value returned. Only if the point was previously + # deleted by the program, restore it; otherwise, leave it + # alone. + + if (fp_equald (wts[i], 0.0D0)) { + if (wts[i] == PDELETE) { + wts[i] = userwts[i] + whydel[i] = NDELETE + } + } else + wts[i] = userwts[i] + } + } + + case HD_K50: + call ic_pstr (ic, "xlabel", "Den + 0.50 * Log (1 - (10 ** -Den))") + xxmin = MIN_DEN + 0.50 * log10 (1.0 - (10. ** (-MIN_DEN))) + xxmax = maxden + 0.50 * log10 (1.0 - (10. ** (-maxden))) + call ic_putr (ic, "xmin", real (xxmin)) + call ic_putr (ic, "xmax", real (xxmax)) + + do i = 1, npts { + dval = Memd[denaf+i-1] + if (dval < 0.0D0 || (fp_equald (dval, 0.0D0))) { + xout[i] = dval + wts[i] = 0.0D0 + whydel[i] = PDELETE + + } else { + xout[i] = dval + 0.50 * log10 (1.0 - (10.** (-dval))) + + # If point had been deleted, find out why. It affects the + # weights value returned. Only if the point was previously + # deleted by the program, restore it; otherwise, leave it + # alone. + + if (fp_equald (wts[i], 0.0D0)) { + if (wts[i] == PDELETE) { + wts[i] = userwts[i] + whydel[i] = NDELETE + } + } else + wts[i] = userwts[i] + } + } + + default: + call eprintf ("Unrecognizable Transform\n") + } + + call sfree (sp) +end diff --git a/noao/imred/dtoi/hdicfit/hdicvshow.x b/noao/imred/dtoi/hdicfit/hdicvshow.x new file mode 100644 index 00000000..e62f6522 --- /dev/null +++ b/noao/imred/dtoi/hdicfit/hdicvshow.x @@ -0,0 +1,155 @@ +include <math/curfit.h> +include "hdicfit.h" + +# IC_VSHOW -- Show fit parameters in verbose mode. + +procedure ic_vshowd (ic, file, cv, x, y, wts, npts, gt) + +pointer ic # ICFIT pointer +char file[ARB] # Output file +pointer cv # Curfit pointer +double x[ARB] # Ordinates +double y[ARB] # Abscissas +double wts[ARB] # Weights +int npts # Number of data points +pointer gt # Graphics tools pointer + +double chisqr, rms +int i, n, deleted, ncoeffs, fd +pointer sp, fit, wts1, coeffs, errors + +int dcvstati(), open() +double ic_rmsd() +errchk open() + +begin + # Do the standard ic_show option, then add on the verbose part. + call ic_show (ic, file, gt) + + if (npts == 0) { + call eprintf ("Incomplete output - no data points for fit\n") + return + } + + # Open the output file. + fd = open (file, APPEND, TEXT_FILE) + + # Determine the number of coefficients and allocate memory. + ncoeffs = dcvstati (cv, CVNCOEFF) + + call smark (sp) + call salloc (coeffs, ncoeffs, TY_DOUBLE) + call salloc (errors, ncoeffs, TY_DOUBLE) + + if (npts == IC_NFIT(ic)) { + # Allocate memory for the fit. + n = npts + call salloc (fit, n, TY_DOUBLE) + call salloc (wts1, n, TY_DOUBLE) + + # Eliminate rejected points and count deleted points. + call amovd (wts, Memd[wts1], n) + if (IC_NREJECT(ic) > 0) { + do i = 1, npts { + if (Memi[IC_REJPTS(ic)+i-1] == YES) + Memd[wts1+i-1] = 0. + } + } + deleted = 0 + do i = 1, n { + if (wts[i] == 0.) + deleted = deleted + 1 + } + + # Get the coefficients and compute the errors. + call dcvvector (cv, x, Memd[fit], n) + call dcvcoeff (cv, Memd[coeffs], ncoeffs) + call dcverrors (cv, y, Memd[wts1], Memd[fit], n, chisqr, + Memd[errors]) + rms = ic_rmsd (x, y, Memd[fit], Memd[wts1], n) + + } else { + # Allocate memory for the fit. + n = IC_NFIT(ic) + call salloc (fit, n, TY_DOUBLE) + call salloc (wts1, n, TY_DOUBLE) + + # Eliminate rejected points and count deleted points. + call amovd (Memd[IC_WTSFIT(ic)], Memd[wts1], n) + if (IC_NREJECT(ic) > 0) { + do i = 1, npts { + if (Memi[IC_REJPTS(ic)+i-1] == YES) + Memd[wts1+i-1] = 0. + } + } + deleted = 0 + do i = 1, n { + if (wts[i] == 0.) + deleted = deleted + 1 + } + + # Get the coefficients and compute the errors. + call dcvvector (cv, Memd[IC_XFIT(ic)], Memd[fit], n) + rms = ic_rmsd (Memd[IC_XFIT(ic)], Memd[IC_YFIT(ic)], + Memd[fit], Memd[wts1], n) + call dcvcoeff (cv, Memd[coeffs], ncoeffs) + call dcverrors (cv, Memd[IC_YFIT(ic)], Memd[wts1], Memd[fit], + n, chisqr, Memd[errors]) + } + + # Print the error analysis. + call fprintf (fd, "total points = %d\n") + call pargi (npts) + call fprintf (fd, "deleted = %d\n") + call pargi (deleted) + call fprintf (fd, "RMS = %7.4g\n") + call pargd (rms) + call fprintf (fd, "square root of reduced chi square = %7.4g\n") + call pargd (sqrt (chisqr)) + + call fprintf (fd, "# \t coefficent\t error\n") + do i = 1, ncoeffs { + call fprintf (fd, "\t%10.4e\t%10.4e\n") + call pargd (Memd[coeffs+i-1]) + call pargd (Memd[errors+i-1]) + } + + # Print x,y pairs and weights + call ic_listxywd (fd, cv, x, y, wts, npts) + + call sfree (sp) + call close (fd) +end + + +# IC_LISTXYW -- List data as x,y pairs on output with their weights. Used +# for verbose show procedure. The untransformed density is also output, +# regardless of what transformation may have been applied. + +procedure ic_listxywd (fd, cv, xvals, yvals, weights, nvalues) + +int fd # File descriptor of output file +pointer cv # Pointer to curfit structure +int nvalues # Number of data values +double xvals[nvalues] # Array of x data values +double yvals[nvalues] # Array of y data values +double weights[nvalues] # Array of weight values + +int i +double dcveval() +include "hdic.com" + +begin + call fprintf (fd,"\n#%15t Density %27t X %39t Yfit %51t LogE %63tWts\n") + + do i = 1, nvalues { + call fprintf (fd, + "%2d %15t%-12.7f%27t%-12.7f%39t%-12.7f%51t%-12.7f%63t%-12.7f\n") + call pargi (i) + call pargd (Memd[den+i-1]) + call pargd (xvals[i]) + call pargd (dcveval (cv, xvals[i])) + call pargd (yvals[i]) + call pargd (weights[i]) + } +end diff --git a/noao/imred/dtoi/hdicfit/mkpkg b/noao/imred/dtoi/hdicfit/mkpkg new file mode 100644 index 00000000..5fc5031f --- /dev/null +++ b/noao/imred/dtoi/hdicfit/mkpkg @@ -0,0 +1,37 @@ +# HDICFIT package. *** Modified from icgfit package for DTOI applications *** + +update: + $update libhdic.a + +libhdic.a: + + hdicclean.x hdicfit.h <pkg/rg.h> + hdicdeviant.x <mach.h> <math/curfit.h> + hdicdosetup.x hdicfit.h <math/curfit.h> + hdicerrors.x hdicfit.h <math/curfit.h> + hdicfit.x hdicfit.h <math/curfit.h> + hdicgaxes.x hdicfit.h <pkg/gtools.h> + hdicgcolon.x hdicfit.h <error.h> <gset.h> + hdicgdelete.x hdicfit.h <gset.h> <mach.h> <pkg/gtools.h> + hdicgfit.x hdicfit.h <error.h> <pkg/gtools.h> <mach.h> <gset.h> + hdicggraph.x hdicfit.h hdic.com <gset.h> <pkg/gtools.h> <mach.h> + hdicgnearest.x <mach.h> <pkg/gtools.h> + hdicgparams.x hdicfit.h <pkg/gtools.h> + hdicgsample.x hdicfit.h <gset.h> <pkg/gtools.h> <pkg/rg.h> + hdicguaxes.x hdic.com hdicfit.h + hdicgundel.x hdicfit.h <gset.h> <mach.h> <pkg/gtools.h> + hdicguser.x + hdicparams.x hdicfit.h + hdicreject.x + hdicshow.x hdicfit.h <pkg/gtools.h> + hdicvshow.x hdicfit.h hdic.com <math/curfit.h> + + hdictrans.x hdicfit.h <mach.h> hdic.com + hdicgvec.x hdicfit.h <mach.h> hdic.com + hdicadd.x hdicfit.h + hdicinit.x hdic.com hdicfit.h <mach.h> + hdicsort.x + userfcn.x <error.h> + hdicgredraw.x <gset.h> + hdicebars.x hdicfit.h hdic.com <pkg/gtools.h> <gset.h> <mach.h> + ; diff --git a/noao/imred/dtoi/hdicfit/userfcn.x b/noao/imred/dtoi/hdicfit/userfcn.x new file mode 100644 index 00000000..d5dba4ed --- /dev/null +++ b/noao/imred/dtoi/hdicfit/userfcn.x @@ -0,0 +1,37 @@ +include <error.h> + +# HD_POWERR -- Construct the basis functions for a power series function. +# Invoked from curfit as a user function. Real version. + +procedure hd_powerr (x, order, k1, k2, basis) + +real x # array of data points +int order # order of polynomial, order = 1, constant +real k1, k2 # normalizing constants - unused +real basis[ARB] # basis functions + +int i + +begin + do i = 1, order + iferr (basis[i] = x ** (i-1)) + call erract (EA_FATAL) +end + + +# HD_POWERD -- Double version of above. + +procedure hd_powerd (x, order, k1, k2, basis) + +double x # array of data points +int order # order of polynomial, order = 1, constant +double k1, k2 # normalizing constants - unused +double basis[ARB] # basis functions + +int i + +begin + do i = 1, order + iferr (basis[i] = x ** (i-1)) + call erract (EA_FATAL) +end |