From 6a9283359bd6d4a6f5a658a9a898d254e75dd71d Mon Sep 17 00:00:00 2001 From: Joern Wilms Date: Mon, 8 Aug 2022 01:25:57 +0200 Subject: [PATCH 1/9] moved tree functions to a new subdirectory --- src/slang/{ => trees}/ntree.sl | 0 src/slang/{ => trees}/ntree_types.sl | 0 2 files changed, 0 insertions(+), 0 deletions(-) rename src/slang/{ => trees}/ntree.sl (100%) rename src/slang/{ => trees}/ntree_types.sl (100%) diff --git a/src/slang/ntree.sl b/src/slang/trees/ntree.sl similarity index 100% rename from src/slang/ntree.sl rename to src/slang/trees/ntree.sl diff --git a/src/slang/ntree_types.sl b/src/slang/trees/ntree_types.sl similarity index 100% rename from src/slang/ntree_types.sl rename to src/slang/trees/ntree_types.sl -- GitLab From 80292609ed4d6aecf936d964b7082829fc631357 Mon Sep 17 00:00:00 2001 From: Jakob Stierhof Date: Thu, 21 Jul 2022 16:24:07 +0200 Subject: [PATCH 2/9] Replace Undefine_Type comparison with __is_initialized test This will take some time probably. New slang version effectively treats all undefined variables as errors. Yes, this renders the Undefine_Type absolutely useless... So, variables can be declared and not initialized, but as soon as the are used things will go south. The proper way to test this is __is_initialized(&potentiall_undefined_variable); Yes, it takes a pointer. Fix here for optimal binning. Oh, and, yes, this is a run time error. That means they will only show up once in use... --- src/data/binning/optimal_binning.sl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/data/binning/optimal_binning.sl b/src/data/binning/optimal_binning.sl index de4ebc5d..e5a0b9f5 100644 --- a/src/data/binning/optimal_binning.sl +++ b/src/data/binning/optimal_binning.sl @@ -186,7 +186,7 @@ private define buildOptimalBinsGroup (id, minCounts) countsId[i] = get_data(fid).value; fakeId[i] = fid; data[i] = did; - if (Undefined_Type == typeof(counts)) + ifnot (__is_initialized(&counts)) counts = Double_Type[length(grid.bin_lo)]; counts += countsId[i]*weights[i]; } -- GitLab From b557a5461796ee19e30e6427bda1d8044b9069a0 Mon Sep 17 00:00:00 2001 From: Jakob Stierhof Date: Thu, 21 Jul 2022 16:56:42 +0200 Subject: [PATCH 3/9] Fixing all Undefine_Type== tests This does not solve the whole issue. Previous version of slang allowed that Undefined_Types can be passed arround as long as they are not actively used. Somthing like this was okay: variable undefined_var; define test (a,b) { return a*100; } test(5, undefined_var); The current devel version of slang throws an error as soon as the undefined type interacts with the stack. So the above code is not valid any more. --- src/FITS/fits_write_tex_table.sl | 2 +- src/Mike/isis_emcee_hammer.sl | 2 +- src/data/lightcurve/lc_from_events.sl | 2 +- src/data/pulsar/arrivaltimes.sl | 6 +++--- src/data/pulsar/arrivaltimes_beta.sl | 6 +++--- src/data/timing/pfold.sl | 2 +- src/fitting/xyfit/array_fit_gauss.sl | 10 +++++----- src/misc/astronomy/orbit_calculator.sl | 2 +- src/misc/sltable.sl | 2 +- src/plot/fplot.sl | 4 ++-- src/slang/arrays/struct_to_assoc_array.sl | 2 +- src/slang/astronomy/jpleph.sl | 6 +++--- src/slang/astronomy/nutation_angles.sl | 2 +- src/slang/astronomy/vsop87.sl | 2 +- src/slang/date/utc2tai.sl | 4 ++-- src/slang/trees/ntree.sl | 2 +- src/slang/trees/ntree_types.sl | 2 +- 17 files changed, 29 insertions(+), 29 deletions(-) diff --git a/src/FITS/fits_write_tex_table.sl b/src/FITS/fits_write_tex_table.sl index 95cd5f0e..70e4169b 100644 --- a/src/FITS/fits_write_tex_table.sl +++ b/src/FITS/fits_write_tex_table.sl @@ -388,7 +388,7 @@ variable inFile; inputFactor=strtrans(inputFactor,"-",""); } inputFactor=sprintf("[$10^{%s}$]",inputFactor); - if(typeof(origName)==Undefined_Type){ + ifnot (__is_initialized(&origName)){ name=oldName+" "+inputFactor+" "+string(units); }else{ name=origName+" "+inputFactor+" "+string(units); diff --git a/src/Mike/isis_emcee_hammer.sl b/src/Mike/isis_emcee_hammer.sl index 9fbe026f..3b2e2c09 100644 --- a/src/Mike/isis_emcee_hammer.sl +++ b/src/Mike/isis_emcee_hammer.sl @@ -184,7 +184,7 @@ public define write_chain( fname, nw, freepar, cdf_scl_lo, cdf_scl_hi, %{{{ variable dlist; list_data(&dlist); - if(typeof(dlist)==Undefined_Type) dlist=" Data not defined at time of file write"; + ifnot(__is_initialized(&dlist)) dlist=" Data not defined at time of file write"; variable hist_comment = ["This file contains the results of a Markov chain analysis", "generated by the emcee hammer subroutine", diff --git a/src/data/lightcurve/lc_from_events.sl b/src/data/lightcurve/lc_from_events.sl index d200c45e..102246ee 100644 --- a/src/data/lightcurve/lc_from_events.sl +++ b/src/data/lightcurve/lc_from_events.sl @@ -56,7 +56,7 @@ define lc_from_events() variable band_min = min(ff.pi); variable band_max = max(ff.pi); - if (typeof(bands) == Undefined_Type) + ifnot (__is_initialized(&bands)) { bands = [band_min, band_max]; } diff --git a/src/data/pulsar/arrivaltimes.sl b/src/data/pulsar/arrivaltimes.sl index a9bfff9b..7d4467be 100644 --- a/src/data/pulsar/arrivaltimes.sl +++ b/src/data/pulsar/arrivaltimes.sl @@ -464,7 +464,7 @@ define atime_merge() { case 1: (a1) = (); } { case 2: (a1,a2) = (); } { help(_function_name()); return; } - ifnot (typeof(a1) == Undefined_Type) a1 = @a1; + if (__is_initialized(&a1)) a1 = @a1; ifnot (a2 == NULL) a2 = @a2; % check on correct parameters @@ -492,8 +492,8 @@ define atime_merge() } else % two structures given { - if (typeof(a1) == Undefined_Type) merg = a2; - else if (typeof(a2) == Undefined_Type) merg = a1; + ifnot (__is_initialized(&a1)) merg = a2; + else ifnot (__is_initialized(&a2)) merg = a1; else { % a2.arrtimes = a2.arrtimes - shift; % correct arrival times diff --git a/src/data/pulsar/arrivaltimes_beta.sl b/src/data/pulsar/arrivaltimes_beta.sl index 87cfc9ed..55181195 100644 --- a/src/data/pulsar/arrivaltimes_beta.sl +++ b/src/data/pulsar/arrivaltimes_beta.sl @@ -485,7 +485,7 @@ define atime_merge_beta() { case 1: (a1) = (); } { case 2: (a1,a2) = (); } { help(_function_name()); return; } - ifnot (typeof(a1) == Undefined_Type) a1 = @a1; + if (__is_initialized(&a1)) a1 = @a1; ifnot (a2 == NULL) a2 = @a2; % check on correct parameters @@ -513,8 +513,8 @@ define atime_merge_beta() } else % two structures given { - if (typeof(a1) == Undefined_Type) merg = a2; - else if (typeof(a2) == Undefined_Type) merg = a1; + ifnot (__is_initialized(&a1)) merg = a2; + else ifnot (__is_initialized(&a2)) merg = a1; else { % a2.arrtimes = a2.arrtimes - shift; % correct arrival times diff --git a/src/data/timing/pfold.sl b/src/data/timing/pfold.sl index 47d46366..bc33de07 100644 --- a/src/data/timing/pfold.sl +++ b/src/data/timing/pfold.sl @@ -367,7 +367,7 @@ define pfold () } {() = (); help(_function_name()); return; } - if (Undefined_Type == typeof(r)) { + ifnot (__is_initialized(&r)) { % event data return periodFoldEvents(t, p;; __qualifiers); } else { diff --git a/src/fitting/xyfit/array_fit_gauss.sl b/src/fitting/xyfit/array_fit_gauss.sl index 9d76de3e..4f5200f7 100644 --- a/src/fitting/xyfit/array_fit_gauss.sl +++ b/src/fitting/xyfit/array_fit_gauss.sl @@ -55,13 +55,13 @@ define array_fit_gauss() ifnot (length(frz)==4) { vmessage("warning (%s): freeze qualifier must be an array of length four",_function_name); frz=[0,0,0,0]; } % starting values if not given variable dx = 0; if (min(x) < 1.) dx = 1.-min(x); % x values must be greater 1 - if (typeof(dy) == Undefined_Type || typeof(dy) == Null_Type) dy = sqrt(y); + if ((not __is_initialized(&dy)) || typeof(dy) == Null_Type) dy = sqrt(y); else ifnot(length(y)==length(dy)) { throw RunTimeError, sprintf("error (%s): y- and dy-arrays must be of equal length",_function_name); } - if (typeof(c) == Undefined_Type || typeof(c) == Null_Type) c = x[where_max(x+dx)]+dx; % center + if ((not __is_initialized(&c)) || typeof(c) == Null_Type) c = x[where_max(x+dx)]+dx; % center else c = c+dx; - if (typeof(s) == Undefined_Type || typeof(s) == Null_Type) s = .5*(max(x)-min(x)); % sigma - if (typeof(a) == Undefined_Type || typeof(a) == Null_Type) a = (max(y)-min(y)); % area - if (typeof(o) == Undefined_Type || typeof(o) == Null_Type) o = 0.; % offset + if ((not __is_initialized(&s)) || typeof(s) == Null_Type) s = .5*(max(x)-min(x)); % sigma + if ((not __is_initialized(&a)) || typeof(a) == Null_Type) a = (max(y)-min(y)); % area + if ((not __is_initialized(&o)) || typeof(o) == Null_Type) o = 0.; % offset % store current fit specifications variable ff = get_fit_fun; variable fm = get_fit_method; diff --git a/src/misc/astronomy/orbit_calculator.sl b/src/misc/astronomy/orbit_calculator.sl index bc8af4af..ffb62032 100644 --- a/src/misc/astronomy/orbit_calculator.sl +++ b/src/misc/astronomy/orbit_calculator.sl @@ -1346,7 +1346,7 @@ define orbit_calculator() %{{{ keys.model = string(model); keys.method = method; keys.eps = string(eps); - if(typeof(seed)!=Undefined_Type) + if(__is_initialized(&seed)) keys = struct_combine( keys, struct{ random_seed = string(seed) } ); if(qualifier_exists("dt")) keys = struct_combine( keys, struct{ dt = string(qualifier("dt")) } ); diff --git a/src/misc/sltable.sl b/src/misc/sltable.sl index c02c4998..0d3b3597 100644 --- a/src/misc/sltable.sl +++ b/src/misc/sltable.sl @@ -395,7 +395,7 @@ private define sltableRound() { limit = struct_field_exists(value,"limit") ? value.limit : 1; nodata = struct_field_exists(value,"nodata") ? value.nodata : 0; value = value.value; - } else if(typeof(mini) != Undefined_Type) { + } else ifnot(__is_initialized(&mini)) { frz = 0; limit = 0; nodata = 0; diff --git a/src/plot/fplot.sl b/src/plot/fplot.sl index 9a9d1930..099f24ab 100644 --- a/src/plot/fplot.sl +++ b/src/plot/fplot.sl @@ -88,7 +88,7 @@ define fplot() { } variable n = fits_read_key(filename,"NUMDAT"); - if (_typeof(dataset) == Undefined_Type) { + ifnot (__is_initialized(&dataset)) { dataset = [1:n]; } %% check if the required datasets are present @@ -100,7 +100,7 @@ define fplot() { } n = length(dataset); - if(_typeof(col) == Undefined_Type) { + ifnot (__is_initialized(&col)) { col = [1:length(dataset)]; } %% check if number of colors match the number of spectra diff --git a/src/slang/arrays/struct_to_assoc_array.sl b/src/slang/arrays/struct_to_assoc_array.sl index 2f0bb12d..39514cfa 100644 --- a/src/slang/arrays/struct_to_assoc_array.sl +++ b/src/slang/arrays/struct_to_assoc_array.sl @@ -34,7 +34,7 @@ define struct_to_assoc_array() % get names of struct fields and initialize array variable fields = get_struct_field_names(s); - if (typeof(type) == Undefined_Type) type = typeof(get_struct_field(s, fields[0])); + ifnot (__is_initialized(&type)) type = typeof(get_struct_field(s, fields[0])); else if (typeof(type) != DataType_Type) { vmessage("error (%s): given type is not of DataType_Type", _function_name); return; } variable a = Assoc_Type[type]; diff --git a/src/slang/astronomy/jpleph.sl b/src/slang/astronomy/jpleph.sl index b01330aa..ce5f90f6 100755 --- a/src/slang/astronomy/jpleph.sl +++ b/src/slang/astronomy/jpleph.sl @@ -710,7 +710,7 @@ define jpl_eph () % (NB: jpl_initeph with no argument defaults to DE430. It % also does all of the error handling % and also search in LHEA_DATA and LHEASOFT) - if (typeof(jpl_ephsave)==Undefined_Type) { + ifnot (__is_initialized(&jpl_ephsave)) { jpl_ephsave=jpl_initeph(); } eph=jpl_ephsave; @@ -1257,7 +1257,7 @@ define planetpos() { case 2: (jdin,objects)=(); } { help(_function_name()); return;} - if (typeof(objects)==Undefined_Type) { + ifnot (__is_initialized(&objects)) { objects=qualifier("object","Sun"); } else { if (qualifier_exists("object")) { @@ -1329,7 +1329,7 @@ define planetpos() % (NB: jpl_initeph with no argument defaults to DE430. It % also does all of the error handling % and also search in LHEA_DATA and LHEASOFT) - if (typeof(jpl_ephsave)==Undefined_Type) { + ifnot (__is_initialized(&jpl_ephsave)) { jpl_ephsave=jpl_initeph(); } eph=jpl_ephsave; diff --git a/src/slang/astronomy/nutation_angles.sl b/src/slang/astronomy/nutation_angles.sl index f0572b9f..45839b35 100755 --- a/src/slang/astronomy/nutation_angles.sl +++ b/src/slang/astronomy/nutation_angles.sl @@ -78,7 +78,7 @@ define nutation_angles() % % read FITS file? % - if (Undefined_Type==typeof(nutation2000a_ls)) { + ifnot (__is_initialized(&nutation2000a_ls)) { variable searchpath=getenv("ISISSCRIPTS_REFPATH"); if (searchpath==NULL and qualifier_exists("iaufile")==0) { throw UsageError,sprintf("%s: iaufile qualifier not given and ISISSCRIPTS_REFPATH not found\n",_function_name()); diff --git a/src/slang/astronomy/vsop87.sl b/src/slang/astronomy/vsop87.sl index 44b22e6e..bab1f6c4 100644 --- a/src/slang/astronomy/vsop87.sl +++ b/src/slang/astronomy/vsop87.sl @@ -197,7 +197,7 @@ variable jdin,planet; } } - if (typeof(vsop87_data)==Undefined_Type) { + ifnot (__is_initialized(&vsop87_data)) { variable searchpath=getenv("ISISSCRIPTS_REFPATH"); if (searchpath==NULL and qualifier_exists("vsopfile")==0) { throw UsageError,sprintf("%s: vsopfile qualifier not given and ISISSCRIPTS_REFPATH not found\n",_function_name()); diff --git a/src/slang/date/utc2tai.sl b/src/slang/date/utc2tai.sl index 8ef94ea3..4c93aa2d 100644 --- a/src/slang/date/utc2tai.sl +++ b/src/slang/date/utc2tai.sl @@ -74,7 +74,7 @@ define utc2tai() help(_function_name()); return; } - if (typeof(tai2utc_data)==Undefined_Type or qualifier_exists("leapfile")) { + if ((not __is_initialized(&tai2utc_data)) or qualifier_exists("leapfile")) { read_leapfile(qualifier("leapfile",getenv("LHEASOFT")+"/refdata/tai-utc.dat")); } @@ -148,7 +148,7 @@ define tai2utc() help(_function_name()); return; } - if (typeof(tai2utc_data)==Undefined_Type or qualifier_exists("leapfile")) { + if ((not __is_initialized(&tai2utc_data)) or qualifier_exists("leapfile")) { read_leapfile(qualifier("leapfile",getenv("LHEASOFT")+"/refdata/tai-utc.dat")); } diff --git a/src/slang/trees/ntree.sl b/src/slang/trees/ntree.sl index bd718656..c1f839ad 100644 --- a/src/slang/trees/ntree.sl +++ b/src/slang/trees/ntree.sl @@ -229,7 +229,7 @@ define ntree() { variable parent = qualifier("parent", NULL); variable index = qualifier("index", NULL); % default structure defining a binary tree - if (typeof(type) == Undefined_Type) + ifnot (__is_initialized(&type)) type = (parent == NULL ? struct { numChilds = 2 } : parent.type); % number of children given else if (typeof(type) == Integer_Type) diff --git a/src/slang/trees/ntree_types.sl b/src/slang/trees/ntree_types.sl index 103ff5cc..99aa3b6e 100644 --- a/src/slang/trees/ntree_types.sl +++ b/src/slang/trees/ntree_types.sl @@ -76,7 +76,7 @@ private define ntree_type_mapQuad_new() { { case 6: (node, c, parent, size, x, y) = (); } { vmessage("error (%s): uncorrect number of arguments", _function_name); return; } - if (parent == NULL && typeof(size) == Undefined_Type) { vmessage("error (%s): size is not specified", _function_name); return; } + if (parent == NULL && (not __is_initialized(&size))) { vmessage("error (%s): size is not specified", _function_name); return; } % extend node structure node = struct_combine(struct { -- GitLab From ef78784461d2fd5785439d574211f8e196120466 Mon Sep 17 00:00:00 2001 From: Ole Koenig Date: Wed, 1 Jun 2022 10:32:23 +0200 Subject: [PATCH 4/9] lag_energy_spectrum: Add coherence, PSD to outfile --- src/data/timing/lag_energy_spectrum.sl | 75 ++++++++++++++++---------- 1 file changed, 47 insertions(+), 28 deletions(-) diff --git a/src/data/timing/lag_energy_spectrum.sl b/src/data/timing/lag_energy_spectrum.sl index 851a94a6..33aae8ac 100644 --- a/src/data/timing/lag_energy_spectrum.sl +++ b/src/data/timing/lag_energy_spectrum.sl @@ -134,8 +134,6 @@ private define write_to_outfile(dat, f_lo, f_hi, keywords, lc_list){ "ENERGY_HI" = qualifier("ehi", Double_Type[N_energies]), "LAG" = dat.lag[*,jj], "LAG_ERR" = dat.errlag[*,jj], - % "LAG_ERR_PROPERR" = dat.errlag_properr[*,jj], - % "LAG_ERR_INGRAM19" = dat.errlag_ingram19[*,jj], "MEAN_IMAG_CPD" = dat.mean_imagcpds[*,jj], "MEAN_IMAG_CPD_ERR" = dat.mean_imagcpd_errs[*,jj], "MEAN_REAL_CPD" = dat.mean_realcpds[*,jj], @@ -146,9 +144,9 @@ private define write_to_outfile(dat, f_lo, f_hi, keywords, lc_list){ }; variable keys = struct_combine(struct{ - F_MIN=f_lo[jj], F_MAX=f_hi[jj], N_SEGMENTS=dat.n_segments[0,jj], - N_FREQUENCIES=dat.n_frequencies[0,jj], F_MIN_DFT=dat.f_min_dft[0,jj], - F_MAX_DFT=dat.f_max_dft[0,jj]}, keywords); + F_MIN = f_lo[jj], F_MAX=f_hi[jj], N_SEGS = dat.n_segments[0,jj], + N_FREQS = dat.n_frequencies[0,jj], FMIN_DFT = dat.f_min_dft[0,jj], + FMAX_DFT = dat.f_max_dft[0,jj]}, keywords); variable hist = struct{ history = [sprintf("Wrote lag-energy spectrum in frequency range %f-%fHz", f_lo[jj], f_hi[jj])], comment = ["Time unit is seconds, energy unit is keV, frequency unit is Hertz", @@ -157,15 +155,33 @@ private define write_to_outfile(dat, f_lo, f_hi, keywords, lc_list){ fits_write_binary_table(fp, extname, sdata, keys, hist); } - %% Write power spectrum into "PSD" extension - variable colnames = String_Type[N_energies+1]; - colnames[0] = "freq"; + %% Write power spectrum and coherence into "COHERENCE" & "PSD" extension + variable psd_colnames = String_Type[N_energies+1]; + variable cof_colnames = String_Type[2*N_energies+1]; + (psd_colnames[0], cof_colnames[0]) = "freq", "freq"; + variable ii; - _for ii (0, N_energies-1, 1) colnames[ii+1] = sprintf("psd%i", ii); - variable psddata = @Struct_Type(colnames); - set_struct_field(psddata, "freq", dat.freq[0]); % frequency grids are all the same (because LCs have same binning) - _for ii (0, N_energies-1, 1) set_struct_field(psddata, colnames[ii+1], dat.psd[ii]); - fits_write_binary_table(fp, "PSD", psddata); + _for ii (0, N_energies-1, 1){ + psd_colnames[ii+1] = sprintf("psd%i", ii); + cof_colnames[ii+1] = sprintf("cof%i", ii); + cof_colnames[N_energies+ii+1] = sprintf("errcof%i", ii); + } + + variable psd_data = @Struct_Type(psd_colnames); + variable cof_data = @Struct_Type(cof_colnames); + + %% Frequency grids are all the same (because LCs have same binning) + set_struct_field(psd_data, "freq", dat.freq[0]); + set_struct_field(cof_data, "freq", dat.freq[0]); + + _for ii (0, N_energies-1, 1){ + set_struct_field(psd_data, psd_colnames[ii+1], dat.psd[ii]); + set_struct_field(cof_data, cof_colnames[ii+1], dat.cof[ii]); + set_struct_field(cof_data, cof_colnames[N_energies+ii+1], dat.errcof[ii]); + } + + fits_write_binary_table(fp, "PSD", psd_data, struct{"N_SEGS"=dat.n_segments[0,0]}); + fits_write_binary_table(fp, "COHERENCE", cof_data, struct{"N_SEGS"=dat.n_segments[0,0]}); fits_close_file(fp); () = system(sprintf("fthedit %s\[0] keyword=N_LES operation=add value=%i comment=\"Number of frequency-resolved lag-energy spectra\"", @@ -244,12 +260,12 @@ define lag_energy_spectrum(lc_list, dimseg) { % % Following keywords are written into each extension: % -% N_SEGMENTS - Number of segments averaged over (foucalc.numavgall) -% N_FREQUENCIES - Number of frequencies averaged over (as in +% N_SEGS - Number of segments averaged over (foucalc.numavgall) +% N_FREQS - Number of frequencies averaged over (as in % "mean of the CPD") % F_MIN/F_MAX - Minimum and maximum frequency as given by % f_lo, f_hi qualifiers -% F_MIN_DFT/F_MAX_DFT - Minimal/Maximal frequency of the frequency-filtered +% FMIN_DFT/FMAX_DFT - Minimal/Maximal frequency of the frequency-filtered % Fourier products - this can be different from % F_MIN/F_MAX due to the discretization of the % *Discrete* Fourier Transform @@ -293,13 +309,15 @@ define lag_energy_spectrum(lc_list, dimseg) { mean_imagcpd_errs = Double_Type[N_energies, N_les], mean_phases = Double_Type[N_energies, N_les], errlag_properr = Double_Type[N_energies, N_les], - % errlag_ingram19 = Double_Type[N_energies, N_les], - freq = Array_Type[N_energies], - psd = Array_Type[N_energies], + f_min_dft = Double_Type[N_energies, N_les], + f_max_dft = Double_Type[N_energies, N_les], n_segments = Integer_Type[N_energies, N_les], n_frequencies = Integer_Type[N_energies, N_les], - f_min_dft = Double_Type[N_energies, N_les], - f_max_dft = Double_Type[N_energies, N_les] + + freq = Array_Type[N_energies], + psd = Array_Type[N_energies], + cof = Array_Type[N_energies], + errcof = Array_Type[N_energies] }; variable lc_ref; @@ -322,9 +340,12 @@ define lag_energy_spectrum(lc_list, dimseg) { rate1 = reference.rate, rate2 = channel_of_interest.rate}, dimseg; normtype=normtype, deadtime=deadtime); + dat.freq[ii] = fouquant.freq; dat.psd[ii] = fouquant.signormpsd2; - + dat.cof[ii] = fouquant.cof12; + dat.errcof[ii] = fouquant.errcof12; + %% Loop through frequency ranges _for jj (0, N_les-1, 1){ if (verbose > 1) vmessage(" ***(%s): Calculating lag-energy spectrum in frequency range %g-%g Hz", @@ -345,26 +366,24 @@ define lag_energy_spectrum(lc_list, dimseg) { %% Coherence and lag calculation on averaged cross-spectrum variable cola = colacal(0.5*(f_lo[jj]+f_hi[jj]), % mid frequency mean(fouquant_f.realcpd12)+mean(fouquant_f.imagcpd12)*1i, - mean(fouquant_f.noicpd12), % not needed for lag+error + mean(fouquant_f.noicpd12), mean(fouquant_f.rawpsd1), mean(fouquant_f.rawpsd2), mean(fouquant_f.noipsd1), mean(fouquant_f.noipsd2), - mean(fouquant_f.sigpsd1), mean(fouquant_f.sigpsd2), % not needed for lag+err + mean(fouquant_f.sigpsd1), mean(fouquant_f.sigpsd2), fouquant_f.numavgall[0]*length(fouquant_f.freq)); % K*M with M=number of frequencies averaged over and K=number of segments (dat.lag[ii,jj], dat.errlag[ii,jj]) = cola.lag, cola.errlag; - + variable m_ReCPD = mean(fouquant_f.realcpd12); variable m_ReCPD_err = mean(fouquant_f.errrealcpd12); variable m_ImCPD = mean(fouquant_f.imagcpd12); variable m_ImCPD_err = mean(fouquant_f.errimagcpd12); %% Calculation from Gaussian error propagation - variable m_phase = atan2(m_ImCPD, m_ReCPD); % radians + variable m_phase = atan2(m_ImCPD, m_ReCPD); % [radians] dat.errlag_properr[ii,jj]=1/(PI*(f_lo[jj]+f_hi[jj]))* sqrt((m_ReCPD/(m_ImCPD^2+m_ReCPD^2))^2*m_ImCPD_err^2 + (m_ImCPD/(m_ImCPD^2+m_ReCPD^2))^2*m_ReCPD_err^2); - %% Error calculation from Ingram 2019, Eq. 19: To be implemented - dat.mean_realcpds[ii,jj] = m_ReCPD; dat.mean_realcpd_errs[ii,jj] = m_ReCPD_err; dat.mean_imagcpds[ii,jj] = m_ImCPD; -- GitLab From e581f19f6262a7191cedf8ca1a8bf5a4473dbfc9 Mon Sep 17 00:00:00 2001 From: Ole Koenig Date: Wed, 1 Jun 2022 10:32:53 +0200 Subject: [PATCH 5/9] segment_lc_for_psd: Add rejected time to out-struct --- src/data/timing/segment_lc_for_psd.sl | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/data/timing/segment_lc_for_psd.sl b/src/data/timing/segment_lc_for_psd.sl index 06ecd1ae..64603b20 100644 --- a/src/data/timing/segment_lc_for_psd.sl +++ b/src/data/timing/segment_lc_for_psd.sl @@ -54,7 +54,8 @@ define segment_lc_for_psd(lc, dt, dimseg) { %{{{ variable segmented_lc = struct{ time=Double_Type[0], rate=Double_Type[0], - error=Double_Type[0] + error=Double_Type[0], + rejected_time, n_segments, total_time }; %% keep a record of the total time and the time thrown away @@ -97,5 +98,9 @@ define segment_lc_for_psd(lc, dt, dimseg) { %{{{ _function_name(), dt*dimseg, dimseg, dt, n_segments); } + segmented_lc.rejected_time = rejected_time; + segmented_lc.total_time = total_time; + segmented_lc.n_segments = n_segments; + return segmented_lc; } %}}} \ No newline at end of file -- GitLab From eb751479f03d18f54216bdc887c1496878522aac Mon Sep 17 00:00:00 2001 From: Joern Wilms Date: Wed, 24 Aug 2022 22:14:29 +0200 Subject: [PATCH 6/9] added various functions for clipping polylines and polygons and two binary tree types --- build/isisscripts_prefix.sl | 2 +- src/slang/geometry/clip_points_polygon.sl | 74 ++ src/slang/geometry/cohen_sutherland.sl | 262 ++++++ src/slang/geometry/greiner_hormann.sl | 452 +++++++++ src/slang/trees/avltree.sl | 1027 +++++++++++++++++++++ 5 files changed, 1816 insertions(+), 1 deletion(-) create mode 100644 src/slang/geometry/clip_points_polygon.sl create mode 100644 src/slang/geometry/cohen_sutherland.sl create mode 100755 src/slang/geometry/greiner_hormann.sl create mode 100644 src/slang/trees/avltree.sl diff --git a/build/isisscripts_prefix.sl b/build/isisscripts_prefix.sl index 81af5628..cc4a676c 100644 --- a/build/isisscripts_prefix.sl +++ b/build/isisscripts_prefix.sl @@ -63,7 +63,7 @@ private define __report_local_paths() define set_local_paths () % This function sets the local paths to find instrument related data and local models -% It can either be done via qualifiers or environemnt variables +% It can either be done via qualifiers or environment variables % The following settings exit: % Data Qualifier Environment % localmodels localmodels XSPEC_LOCAL_MODELS,XSPEC_TABLE_MODELS % semi-colon seperated list of paths diff --git a/src/slang/geometry/clip_points_polygon.sl b/src/slang/geometry/clip_points_polygon.sl new file mode 100644 index 00000000..92b3caa3 --- /dev/null +++ b/src/slang/geometry/clip_points_polygon.sl @@ -0,0 +1,74 @@ +define clip_points_polygon() +%!%+ +%\function{clip_points_polygon} +%\synopsis{Clip points against a polygon} +%\usage{(xc,yc)=clip_points_polygon(x,y,xp,yp)} +%\altusage{clipped=cohen_sutherland{points,poly}} +%\qualifiers{ +% \qualifier{evenodd}{use the even-odd method to determine } +% \qualifier{crossing}{use the crossing number method} +% \qualifier{winding}{use the winding number method (the default)} +%} +% +%\description +%This function clips points defined by (x,y), where x and y can be arrays, +%against a closed polygon defined by the arrays (xp, yp) and returns the +%clipped points (xc,yc) where xc, yc are arrays (which can be empty +%if all points are outside of the polygon). Here, a closed polygon +%means that xp[0]==xp[-1] and yp[0]==yp[-1]. +% +%Alternatively, the points and polygon can be defined as structs, +%where points=struct{ x=[], y=[] } and where the polygon is defined +%as poly=struct{x=xp,y=yp} +% +%The qualifiers define what to consider the "inside" of the polygon. +% +%\seealso{clip_points_rectangle, greiner_hormann, point_in_polygon} +% +%!%- +{ + variable src,clp; + variable xx,yy,vx,vy; + + variable stru=0; + if (_NARGS == 2 ) { + (src,clp)=(); + stru=1; + xx=src.x; + yy=src.y; + vx=clp.x; + vy=clp.y; + } else { + if (_NARGS == 4 ) { + (xx,yy,vx,vy)=(); + } else { + throw UsageError,"clip_points_polygon: Use either 2 or 4 arguments\n"; + } + } + + variable rx={}; + variable ry={}; + variable i; + _for i(0,length(xx)-1,1) { + if (point_in_polygon(xx[i],yy[i],vx,vy;;__qualifiers())) { + list_append(rx,xx[i]); + list_append(ry,yy[i]); + } + } + + if (length(rx)==0) { + if (stru) { + return struct {x=Double_Type[0],y=Double_Type[0]}; + } + return (Double_Type[0],Double_Type[0]); + } + + list_to_array(xx); + list_to_array(yy); + + if (stru) { + return struct{x=xx,y=yy}; + } + + return (xx,yy); +} diff --git a/src/slang/geometry/cohen_sutherland.sl b/src/slang/geometry/cohen_sutherland.sl new file mode 100644 index 00000000..4cc0b75d --- /dev/null +++ b/src/slang/geometry/cohen_sutherland.sl @@ -0,0 +1,262 @@ + +private variable CS_inside=0; % 0000 +private variable CS_left=1; % 0001 +private variable CS_right=2; % 0010 +private variable CS_bottom=4; % 0100 +private variable CS_top=8; % 1000 + +private define cs_outcode(x,y,xmin,ymin,xmax,ymax) { + % + % return location of x,y against box xmin,ymin,xmax,ymax + variable code=CS_inside; + if (xxmax) { + code |=CS_right; + } + } + if (yymax) { + code |=CS_top; + } + } + return code; +} + + +define cohen_sutherland() +%!%+ +%\function{cohen_sutherland} +%\synopsis{Clip a line against a rectangle} +%\usage{(xc0,yc0,xc1,yc1)=cohen_sutherland(x0,y0,x1,y1,xmin,ymin,xmax,ymax)} +%\altusage{clipped=cohen_sutherland(line,box)} +%\description +%This function clips a line defined by the points (x0,y0) and (x1,y1) against +%a box defined by the corner points (xmin,ymin) and (xmax,ymax) and returns +%the clipped line (xc0,yc0) -- (xc1,yc1). The clipped coordinates are set +%to _NaN if the line misses the box. +% +%The line segment and clipping box can be either defined directly by giving the +%coordinates or as structs, In the latter case, the line segment +%is defined as line= struct{ x=[x0,x1], y=[y0,y1] } and the +%rectangular box is defined either as box=struct{x=[xmin,xmax],y=[ymin,ymax]} +%or as box=struct{xmin=xmin,ymin=ymin,xmax=xmax,ymax=ymax}. +% +%\seealso{clip_points_rectangle,clip_points_polygon, clip_polyline_rectangle,greiner_hormann} +% +%!%- +{ + variable x0,y0,x1,y1,xmin,ymin,xmax,ymax; + + variable struret=0; + + if ( _NARGS == 2 ) { + variable p,box; + (p,box)=(); + x0=p.x[0]; y0=p.y[0]; + x1=p.x[1]; y1=p.y[1]; + + if (struct_field_exists(box,"x")) { + xmin=box.x[0]; xmax=box.x[1]; + ymin=box.y[0]; ymax=box.y[1]; + } else { + % shp convention for bounding boxes + xmin=box.xmin; ymin=box.ymin; + xmax=box.xmax; ymax=box.ymax; + } + + struret=1; + + } else { + if (_NARGS==8) { + (x0,y0,x1,y1,xmin,ymin,xmax,ymax)=(); + } else { + throw UsageError,"cohen_sutherland: Use either 2 or 8 arguments\n"; + } + } + + variable outcode0=cs_outcode(x0,y0,xmin,ymin,xmax,ymax); + variable outcode1=cs_outcode(x1,y1,xmin,ymin,xmax,ymax); + + variable accept=0; + + while(1) { + if ( not (outcode0 | outcode1)) { + % bitwise or is 0: both points are inside the window + accept=1; + break; + } else { + % bitwise and is not 0: both points are outside the window + if (outcode0 & outcode1) { + break; + } else { + % both tests failed: calculate the line segment + % to clip from an outside point to an intersection + % with the clip edge; + + variable x,y; + variable outcodeOut = outcode1 > outcode0 ? outcode1 : outcode0; + + % find the intersection point + if (outcodeOut & CS_top) { %point is above clip window + x=x0+(x1-x0)*(ymax-y0)/(y1-y0); + y=ymax; + } else { + if (outcodeOut & CS_bottom) { % point is below clip window + x=x0+(x1-x0)*(ymin-y0)/(y1-y0); + y=ymin; + } else { + if (outcodeOut & CS_right) { % right of clip window + x=xmax; + y=y0+(y1-y0)*(xmax-x0)/(x1-x0); + } else { + if (outcodeOut & CS_left) { % left of clip window + x=xmin; + y=y0+(y1-y0)*(xmin-x0)/(x1-x0); + } + } + } + } + + % move outside point to intersection point + % and move to next pass + if (outcodeOut==outcode0) { + x0=x; + y0=y; + outcode0=cs_outcode(x0,y0,xmin,ymin,xmax,ymax); + } else { + x1=x; + y1=y; + outcode1=cs_outcode(x1,y1,xmin,ymin,xmax,ymax); + } + } + } + } + + if (accept) { + if (struret) { + return struct{x=[x0,x1],y=[y0,y1]}; + } + return (x0,y0,x1,y1); + } + + % segment is not in the clip window + if (struret) { + return struct{x=Double_Type[0],y=Double_Type[0]}; + } + return (_NaN,_NaN,_NaN,_NaN); + +} + +define clip_polyline_rectangle(src,box) { +%!%+ +%\function{clip_polyline_rectangle} +%\synopsis{Clip a polyline against a rectangle} +%\usage{clipped=clip_polyline_rectangle(poly,box)} +%\description +%This function clips the polyline poly=struct{x=[],y=[]}, where x and y are arrays containing +%the points of the polyline, against a rectangle defined by the corner points (xmin,ymin) and +%(xmax,ymax). It returns a list of structs{x,y}, where each list element contains the points of +%a segment of the polygon that is inside the box. +% +%\seealso{clip_points_polygon, clip_polyline_rectangle, cohen_sutherland, greiner_hormann} +% +%!%- + + variable i; + + variable polyli={}; + variable xx={}; + variable yy={}; + + _for i(0,length(src.x)-2,1) { + variable seg=cohen_sutherland(struct{x=[src.x[i],src.x[i+1]], + y=[src.y[i],src.y[i+1]]},box); + % is the segment within the box? + if (length(seg.x)>0) { + % append? + if (length(xx)>0) { + % start of new segment different from old end point? + if ((xx[-1]!=seg.x[0]) || yy[-1]!=seg.y[0]) { + xx=list_to_array(xx); + yy=list_to_array(yy); + list_append(polyli,struct{x=xx,y=yy}); + xx={}; + yy={}; + } + } + list_append(xx,seg.x[0]); list_append(xx,seg.x[1]); + list_append(yy,seg.y[0]); list_append(yy,seg.y[1]); + } + } + + if (length(xx)>0) { + xx=list_to_array(xx); + yy=list_to_array(yy); + list_append(polyli,struct{x=xx,y=yy}); + } + + return polyli; +} + +define clip_points_rectangle() +%!%+ +%\function{clip_points_rectangle} +%\synopsis{Clip points against a rectangle} +%\usage{(xc,yc)=clip_points_rectangle(x,y,xmin,ymin,xmax,ymax)} +%\altusage{clipped=clip_points_rectangle(points,box)} +%\description +%This function clips the points defined by (x,y), where x and y can be arrays, +%against a rectangle defined by the corner points (xmin,ymin) and (xmax,ymax) and +%returns the clipped points (xc,yc) where xc, yc are arrays (which can be empty +%if all points are outside of the rectangle). +% +%Alternatively, the points and clipping rectangle can be defined as structs, +%where points=struct{ x=[], y=[] } and where the clipping rectangle is defined +%either as box=struct{x=[xmin,xmax],y=[ymin,ymax]} +%or as box=struct{xmin=xmin,ymin=ymin,xmax=xmax,ymax=ymax}. +% +% +%\seealso{clip_points_polygon, clip_polyline_rectangle, cohen_sutherland, greiner_hormann} +% +%!%- +{ + variable xmin,xmax,ymin,ymax; + variable struret=0; + variable xi,yi; + + if ( _NARGS == 2 ) { + variable src,box; + (src,box)=(); + xi=src.x; + yi=src.y; + + if (struct_field_exists(box,"x")) { + xmin=box.x[0]; xmax=box.x[1]; + ymin=box.y[0]; ymax=box.y[1]; + } else { + % shp convention for bounding boxes + xmin=box.xmin; ymin=box.ymin; + xmax=box.xmax; ymax=box.ymax; + } + + struret=1; + } else { + if (_NARGS==6) { + (xi,yi,xmin,ymin,xmax,ymax)=(); + } else { + throw UsageError,"clip_points: Use either 2 or 6 arguments\n"; + } + } + + variable ndx=where(xmin<=xi0.); + + variable S02x=P0.x-P2.x; + variable S02y=P0.y-P2.y; + + variable snumer=S10x*S02y-S10y*S02x; + + if ((snumer<0)==denomPositive) { + return 0; % no intersection + } + + variable tnumer=S32x*S02y-S32y*S02x; + if ((tnumer<0.)==denomPositive) { + return 0; % no intersection + } + + if (((snumer>denom)==denomPositive) or ((tnumer>denom)==denomPositive)) { + return 0; % no intersection + } + + @alphaP=tnumer/denom; + @alphaQ=snumer/denom; + return 1; +} + +private define new_greiner_vertex() { + % + % private function: define a new vertex point for the + % Greiner and Hormann algorithm + % + variable x,y,id,alpha,inter; + inter=0; + alpha=0.; + if (_NARGS==3) { + (x,y,id)=(); + } + if (_NARGS==4) { + (x,y,id,alpha)=(); + inter=1; + } + + return struct { + x=x, % x-coordinate + y=y, % y-coordinate + id=id, % identifier + next=NULL, % pointer to next vertex point + prev=NULL, % pointer to current vertex point + intersect=inter, % 0 not an intersect, 1 intersect + neighbor=NULL, % pointer to this point in other array + alpha=alpha, % storage for distance from last non-intersection point + entry=-1, % 1 if entry point, 0 if not, -1 undefined + processed=0 % true if vertex has been taken care of + }; +} + +private define new_greiner_vertexlist(P) { + % build up a vertex list for the Greiner-Hormann algorithm + variable i; + variable sstart=new_greiner_vertex(P.x[0],P.y[0], P.id[0]); + variable ssold=sstart; + _for i(1,length(P.x)-1,1) { + variable ssnew=new_greiner_vertex(P.x[i],P.y[i], P.id[i]); + ssnew.prev=ssold; + ssold.next=ssnew; + ssold=ssnew; + } + sstart.prev=ssold; + ssnew.next=sstart; + return sstart; +} + +private define greiner_insert_between(P,Q,ins) { + % insert vertex np between non-intersecting points P and Q + % nb: does not check that P, Q are non intersection points + % for speed reasons + + % np must be an intersection vertex with alpha set + + % find correct insertion segment + variable pp=P; + while (pp.next != Q && pp.next.alpha1) { + throw UsageError,sprintf("Usage error in %s: only one qualifier of union, intersection, and without can be set!\n", _function_name()); + } + + % default: clip S with C (i.e., get the intersection of the polygons) + if (sm==0) { + inter=1; + } + + % build up vertex list of the subject and the clip polygons + variable S,C; + if (_NARGS==2) { + variable Stmp,Ctmp; + (Stmp,Ctmp)=(); + S=struct { x=@Stmp.x, y=@Stmp.y, id=Char_Type[length(Stmp.x)]}; + C=struct { x=@Ctmp.x, y=@Ctmp.y, id=1+Char_Type[length(Ctmp.x)]}; + } else { + variable sx,sy,cx,cy; + (sx,sy,cx,cy)=(); + S=struct { x=sx[*], y=sy[*], id=Char_Type[length(sx)] }; + C=struct { x=cx[*], y=cy[*], id=1+Char_Type[length(cx)] }; + } + + if (S.x[0]!=S.x[-1] || S.y[0]!=S.y[-1] ) { + throw UsageError,sprintf("Usage error in %s: Source polygon is not closed!\n", _function_name()); + } + + if (C.x[0]!=C.x[-1] || C.y[0]!=C.y[-1] ) { + throw UsageError,sprintf("Usage error in %s: Clip polygon is not closed!\n", _function_name()); + } + + if (length(S.x)!=length(S.y)) { + throw UsageError,sprintf("Usage error in %s: Source x and y arrays must be of equal length\n",_function_name()); + } + + if (length(C.x)!=length(C.y)) { + throw UsageError,sprintf("Usage error in %s: Clip x and y arrays must be of equal length\n",_function_name()); + } + + if (qualifier_exists("perturb")) { + S.x+=1e-8*urand(length(S.x)); + S.y+=1e-8*urand(length(S.x)); + + S.x[-1]=S.x[0]; + S.y[-1]=S.y[0]; + } + + + variable src=new_greiner_vertexlist(S); + variable clp=new_greiner_vertexlist(C); + + % we'll be returning a list of polygons + variable polylist={}; + + % + % Step 1: find all intersection points between the polygons and + % insert them into the polygons + % + variable numinter=0; + variable ss=src; + do { + % find next vortex in src that is not an intersection + variable snxt=ss; + do { + snxt=snxt.next; + } while (snxt.intersect==1); + + variable cc=clp; + do { + % find next vortex that is not an intersection + variable cnxt=cc; + do { + cnxt=cnxt.next; + } while (cnxt.intersect==1); + + % Is there an intersection? + variable alphaP,alphaQ; + if (line_intersect(ss,snxt,cc,cnxt,&alphaP,&alphaQ)==1) { + numinter++; + % insert new intersection point into both polygons + variable newx=ss.x+alphaP*(snxt.x-ss.x); + variable newy=ss.y+alphaP*(snxt.y-ss.y); + + variable np1=new_greiner_vertex(newx,newy,-1,alphaP); + variable np2=new_greiner_vertex(newx,newy,-1,alphaQ); + np1.neighbor=np2; + np2.neighbor=np1; + greiner_insert_between(ss,snxt,np1); + greiner_insert_between(cc,cnxt,np2); + } + cc=cnxt; + } while (cc!=clp); + ss=snxt; + + } while (ss!=src); + + % + % Deal with the case that no intersections were found + % + if (numinter==0) { + % is S inside C or outside? + variable insi=point_in_polygon(S.x[0],S.y[0],C.x,C.y); + + if (insi) { + % S is inside C + if (uni) { + list_append(polylist,struct{x=C.x[*],y=C.y[*],id=C.id[*]}); + } + + % S without C gives empty list + if (inter) { + list_append(polylist,struct{x=S.x[*],y=S.y[*],id=C.id[*]}); + } + } else { + % S is outside C + + % union and without C + if (uni or without) { + list_append(polylist,struct{x=S.x[*],y=S.y[*],id=S.id[*]}); + } + + % intersection gives empty list + } + return polylist; + } + + % + % Phase 2: identify the entry and exit points; + % depending on the final operation desired, the + % initialization is done in a slightly different way + % + + % these different starting values will result in our three + % cases + variable sentry,centry; + if (uni) { + sentry=0; + centry=0; + } + if (inter) { + sentry=1; + centry=1; + } + if (without) { + sentry=0; + centry=1; + } + + % is starting point for ss inside of C? + ss=src; + variable entry=sentry xor point_in_polygon(ss.x,ss.y,C.x,C.y); + while (ss.next!=src) { + if (ss.intersect) { + ss.entry=entry; + entry=not entry; + } + ss=ss.next; + } + + % is starting point for cc inside of S? + cc=clp; + entry=centry xor point_in_polygon(cc.x,cc.y,S.x,S.y); + while (cc.next!=clp) { + if (cc.intersect) { + cc.entry=entry; + entry=not entry; + } + cc=cc.next; + } + + % + % Step 3: Go through polygons and build up the new structure + % + variable unprocessed=1; + do { + + % + % go to the first unchecked intersection in the list + % + variable xx={}; + variable yy={}; + variable id={}; + + cc=src; + while ( (cc.intersect==0 || cc.processed==1) && cc.next!=src) { + cc=cc.next; + } + do { + cc.processed=1; + if (cc.entry==1) { + do { + list_append(xx,cc.x); + list_append(yy,cc.y); + list_append(id,cc.id); + cc.processed=1; + cc=cc.next; + } while (cc.intersect==0); + } else { + do { + list_append(xx,cc.x); + list_append(yy,cc.y); + list_append(id,cc.id); + cc.processed=1; + cc=cc.prev; + } while (cc.intersect==0); + } + cc.processed=1; + cc=cc.neighbor; + } while (cc.processed==0); + + % save the current polygon on the output list + xx=list_to_array(xx); + yy=list_to_array(yy); + id=list_to_array(id); + list_append(polylist,struct{x=[xx,xx[0]],y=[yy,yy[0]],id=[id,id[0]]}); + + cc=src; + unprocessed=0; + while (cc.next!=src) { + if (cc.intersect==1 && cc.processed==0) { + unprocessed++; + } + cc=cc.next; + } + + % is there still an unprocessed intersection in the list? + cc=src; + unprocessed=0; + while (unprocessed==0 && cc.next!=src) { + if (cc.intersect==1 && cc.processed==0) { + unprocessed=1; + } + cc=cc.next; + } + } while (unprocessed==1); + + return polylist; +} diff --git a/src/slang/trees/avltree.sl b/src/slang/trees/avltree.sl new file mode 100644 index 00000000..8f4af529 --- /dev/null +++ b/src/slang/trees/avltree.sl @@ -0,0 +1,1027 @@ +% +% Implementation of binary search and AVL trees +% + +% Basic type definitions +private variable TREE_LEFT=0; +private variable TREE_RIGHT=1; + +% +% basic key compare function +% (can be overridden if the tree contains more complex data structures) +% +private define tree_key_compare(k1,k2) { + return (k1 < k2) ? -1 : ( (k1==k2) ? 0 : +1); +} + +% +% Binary tree specific functions and structures +% +private define new_BinarySearchTreeNode() { + variable T; + variable data=NULL; + if (_NARGS==1) { + T=(); + } + if (_NARGS==2) { + (T,data)=(); + } + + return struct { % node of a binary search tree + child={NULL,NULL}, % value is NULL if child is empty + data=data + }; +} + +private define BinarySearchTreeInsert_helper(T,key) { + % insert a node containing key into the binary tree T + % this could also be programmed recursively, but iteratively is + % faster. This returns the newly inserted key. + + % use new_node function such that this will also work on + % binary search trees derived from the BinarySearchTree + variable N=T.new_node(key); + + % empty tree: insert at root and be done with it + if (typeof(T.root)!=Struct_Type) { + T.root=N; + return N; + } + + % existing tree - iterative traversal to correct position + variable no=T.root; + while (1) { + variable cmp=T.key_compare(key,no.data); + % key < no.data ? + if (cmp<0) { + % empty child? - if yes, then insert + if (no.child[TREE_LEFT]==NULL) { + no.child[TREE_LEFT]=N; + return N; + } + no=no.child[TREE_LEFT]; + } else { + % key > no.data ? + if (cmp>0) { + % empty child? - if yes, then insert + if (no.child[TREE_RIGHT]==NULL) { + no.child[TREE_RIGHT]=N; + return N; + } + no=no.child[TREE_RIGHT]; + } else { + % key already exists, throw an error + % FIXME: add graceful exit option + throw UsageError,sprintf("insert: Key %S already exists!\n",key); + } + } + } + + % should never reach this point + throw UsageError,sprintf("This should never happen\n"); +} + +private define BinarySearchTreeInsert(T,key) { + % insert a node containing key into the binary tree T + ()=BinarySearchTreeInsert_helper(T,key); +} + +private define BinarySearchTreeMinimum(N) { + % return node with the smallest key in a binary search tree + % if qualifier "node" is set returns the node, otherwise the key + % + % N can be a node or a tree + + % case of tree + if (struct_field_exists(N,"root")) { + N=N.root; + } + + while (N.child[TREE_LEFT]!=NULL) { + N=N.child[TREE_LEFT]; + } + if (qualifier_exists("node")) { + return N; + } + return N.data; +} + +define BinarySearchTreeMaximum(N) { + % return the largest key in a binary search tree + % if qualifier "node" is set returns the node, otherwise the key + + % case of tree + if (struct_field_exists(N,"root")) { + N=N.root; + } + + while (N.child[TREE_RIGHT]!=NULL) { + N=N.child[TREE_RIGHT]; + } + if (qualifier_exists("node")) { + return N; + } + return N.data; +} + +% iterative deletion of a key in a binary search tree +define BinarySearchTreeDeleteKey_helper_iterative(T,key,N); +define BinarySearchTreeDeleteKey_helper_iterative(T,key,N) { + % + % deletion of node with key K from binary tree (version 1) + % + if (N==NULL) { + return NULL; + } + + % + % traverse tree + % + variable cmp=T.key_compare(key,N.data); + variable parent; + while (N!=NULL && cmp!=0) { + parent=N; + if (cmp<0) { + N=N.child[TREE_LEFT]; + } else { + N=N.child[TREE_RIGHT]; + } + + % not found + if (N==NULL) { + return N; + } + cmp=T.key_compare(key,N.data); + } + + % N is now the key to be deleted + + if (N.child[TREE_LEFT]==NULL) { + if (N.child[TREE_RIGHT]==NULL) { + N=NULL; % kill data + return N; + } + N=N.child[TREE_RIGHT]; + return N; + } + + if (N.child[TREE_RIGHT]==NULL) { + N=N.child[TREE_LEFT]; + return N; + } + + % node has two children: replace current node with successor data + % and delete successor + + % find successor node + variable succ=BinarySearchTreeMinimum(N.child[TREE_RIGHT];node); + % copy successor data + N.data=succ.data; + + % delete inorder successor + N.child[TREE_RIGHT]=BinarySearchTreeDeleteKey_helper_iterative(T,N.data,N.child[TREE_RIGHT]); + + return N; +} + +% +% full recursive tree deletion (needed for the AVL tree) +% +define BinarySearchTreeDeleteKey_helper_recursive(T,key,N); + +define BinarySearchTreeDeleteKey_helper_recursive(T,key,N) { + % no node at current position, end recursion + if (N==NULL) { + return NULL; + } + + % recursive deletion + variable cmp=T.key_compare(key,N.data); + if (cmp<0) { + % key < N.data + N.child[TREE_LEFT]=BinarySearchTreeDeleteKey_helper_recursive(T,key,N.child[TREE_LEFT]); + } else { + if (cmp>0) { + % key > N.data + N.child[TREE_RIGHT]=BinarySearchTreeDeleteKey_helper_recursive(T,key,N.child[TREE_RIGHT]); + } else { + % all ifs below: found it! + % + + if (N.child[TREE_LEFT]==NULL && N.child[TREE_RIGHT]==NULL) { + % no child present, just delete the thing + N=NULL; + } else { + if (N.child[TREE_LEFT]==NULL) { + N=N.child[TREE_RIGHT]; + } else { + if (N.child[TREE_RIGHT]==NULL) { + N=N.child[TREE_LEFT]; + } else { + % delete node with two children + + % find inorder successor of current node + variable succ=BinarySearchTreeMinimum(N.child[TREE_RIGHT];node); + % save the data of the successor + N.data=succ.data; + % delete successor + N.child[TREE_RIGHT]=BinarySearchTreeDeleteKey_helper_recursive(T,N.data,N.child[TREE_RIGHT]); + } + } + } + } + } + return N; +} + +define BinarySearchTreeDeleteKey(T,key) { + ()=BinarySearchTreeDeleteKey_helper_iterative(T,key,T.root); +} + +define BinarySearchTreeDeleteNode(T,N) { + ()=BinarySearchTreeDeleteKey_helper_iterative(T,N.key,T.root); +} + +% find node for which data==key +% returns the Node of Null, if the key does not exist +define BinarySearchTreeSearch(T,key) { + variable N=T.root; + while (N!=NULL && N.data != key ) { + if (key < N.data) { + N=N.child[TREE_LEFT]; + } else { + N=N.child[TREE_RIGHT]; + } + } + + return N; +} + +% return true if key is contained in T +private define BinarySearchTreeKeyExists(T,key) { + variable N=BinarySearchTreeSearch(T,key); + return (N!=NULL); +} + +private define BinarySearchTreeNodeSuccessor(T,N) +% find successor node of N in tree T +% (T isn't really needed here, but this helps the tree object) +{ + if (N.child[TREE_RIGHT]!=NULL) { + return BinarySearchTreeMinimum(N.child[TREE_RIGHT];node); + } + variable y=N.parent; + while (y!=NULL && N==y.child[TREE_RIGHT]) { + N=y; + y=y.parent; + } + return y; +} + +private define BinarySearchTreeNodePredecessor(T,N) +% find predecessor node of N in tree T +% (T isn't really needed here, but this helps the tree object) +{ + if (N.child[TREE_LEFT]!=NULL) { + return BinarySearchTreeMinimum(N.child[TREE_LEFT];node); + } + variable y=N.parent; + while (y!=NULL && N==y.child[TREE_LEFT]) { + N=y; + y=y.parent; + } + return y; +} + +private define BinarySearchPrintNode(N) { + ()=printf("%S\n",N.data); +} + +private define BinarySearchTraversal(); + +private define BinarySearchTraversal() { + variable N,func; + + if (_NARGS==1) { + N=(); + func=&BinarySearchPrintNode; + } + + if (_NARGS==2) { + (N,func)=(); + } + + if (N==NULL) { + return; + } + + if (struct_field_exists(N,"root")) { + N=N.root; + if (N==NULL) { + return; + } + } + + if (qualifier_exists("preorder")) { + (@func)(N;;__qualifiers()); + BinarySearchTraversal(N.child[TREE_LEFT],func;;__qualifiers()); + BinarySearchTraversal(N.child[TREE_RIGHT],func;;__qualifiers()); + return; + } + + if (qualifier_exists("postorder")) { + BinarySearchTraversal(N.child[TREE_LEFT],func;;__qualifiers()); + BinarySearchTraversal(N.child[TREE_RIGHT],func;;__qualifiers()); + (@func)(N;;__qualifiers()); + return; + } + + % default is inorder traversal + BinarySearchTraversal(N.child[TREE_LEFT],func;;__qualifiers()); + (@func)(N;;__qualifiers()); + BinarySearchTraversal(N.child[TREE_RIGHT],func;;__qualifiers()); + return; + +} + +private define BinarySearchTreeNodePrint(); + +private define BinarySearchTreeNodePrint() { + variable N,sp; + + (N,sp)=(); + + if (N==NULL) { + return; + } + + BinarySearchTreeNodePrint(N.child[TREE_RIGHT],sp+" "); + if (struct_field_exists(N,"height")) { + ()=printf("%s%S-%i\n",sp,N.data,N.height); + } else { + ()=printf("%s%S\n",sp,N.data); + } + BinarySearchTreeNodePrint(N.child[TREE_LEFT],sp+" "); +} + +private define BinarySearchTreePrint(T) { + BinarySearchTreeNodePrint(T.root,""); +} + +variable BinarySearchTree_Struct=struct { + root=NULL, + insert=&BinarySearchTreeInsert, + print=&BinarySearchTreePrint, + new_node=&new_BinarySearchTreeNode, + key_compare=&tree_key_compare, + search=&BinarySearchTreeSearch, + delete=&BinarySearchTreeDeleteKey, + delete_node=&BinarySearchTreeDeleteNode, + exists=&BinarySearchTreeKeyExists, + min=&BinarySearchTreeMinimum, + max=&BinarySearchTreeMaximum, + successor=&BinarySearchTreeNodeSuccessor, + predecessor=&BinarySearchTreeNodePredecessor, + traversal=&BinarySearchTraversal +}; + +define BinarySearchTree() +%!%+ +%\function{BinarySearchTree} +%\synopsis{Creates a new binary search tree} +%\usage{Struct_Type=BinarySearchTree();} +%\altusage{Struct_Type=BinarySearchTree(cmpfunction);} +%\description +%A binary search tree is an object that permits fast +%searches and ordered operations on ordinal data types (i.e., +%data that can be sorted). +% +%The tree consists of nodes which are structs containing +%the data (or "key") in tag "data" and references to nodes with +%keys with values less than data in child[0] and to keys with +%values larger than data in child[1]. In general, users of the +%data structure should not worry about the structure of the +%nodes and use the functions described below for all tree +%operations. +% +%Searches in a binary search tree can be fast (O(log N)) if the +%search tree is balanced, that is, if there is a similar number +%of nodes below each node (if that makes sense...). This is typically +%the case if random (unsorted) data are inserted into the tree, and +%very much not the case if ordered data are inserted. In this case +%the search degenerates to O(N). The isisscripts provide a special +%binary search tree called AVLTree that has the same accessor functions +%as BinarySearchTree but ensures that the tree is height-balanced (at some +%small additional cost for key insertions and deletions). +% +%The BinarySearchTree and the AVLTree provide the following functions: +% insert(key) - insert a data element key into the binary search tree +% delete(key) - remove a data element key from the binary search tree +% exists(key) - return true if the search tree contains key +% search(key) - return the node containing the key (or NULL if +% they key does not exist) +% min - return the minimum key of the search tree (or a subtree), or +% the node containing the minimum key if the node qualifier +% is set. +% max - return the maximum key of the search tree (or a subtree), or +% the node containing the maximum key if the node qualifier +% is set. +% successor - return the successor node to a node (i.e., the node +% containing the next largest data value in the tree). +% predecessor - return the predecessor node to a node (i.e., the node +% containing the next smaller data value in the tree). +% print - print the structure of the tree +% traversal(&func;quals) - traverse the tree, i.e., execute func for all +% keys of the tree. The function func is called with each +% node N (i.e., the data are found in N.data) and all +% qualifiers given to traversal. The order in which the +% nodes are traversed is defined by the following qualifiers: +% preorder: operate on the node, then on the left children, +% then on the right children +% postorder: operate first on the left children, then on the +% right children, then on the node +% inorder (the default): operate on the left children, then +% on the node, then on the right children +% Note that it is permitted that func modifies the tree (e.g., +% insert new elements, delete tree elements and so on; it is +% best practice to use a qualifier that specifies the tree to +% let func know about it). Depending on the traversal order +% the modified elements may or may not be operated on as part +% of the traversal. +% +% +%The default setup of the binary search tree works on all data types +%where the operators <, ==, and > have been defined. For other data +%types, initialize the tree with a reference to a comparison function +%of the style define key_compare(k1,k2) where k1 and k2 are of the +%data type contained in the search tree and which returns a negative +%number if k1k2 (this +%is the same definition as the one used by s-lang's sort function or +%by the strcmp function). +% +%\example +% % sort some numbers +% % (obviously a better way would be array_sort...) +% +% define arrapp(N) { +% % append key of N to list +% variable list=qualifier("list"); +% list_append(list,N.data); +% } +% +% variable arr=[1,4,6,3,5,7,8,2,9]; +% variable t=BinarySearchTree(); +% variable i; +% foreach i (arr) { +% t.insert(i); +% } +% +% % print the tree +% t.print(); +% +% % delete element with key 8 +% t.delete(8); +% +% variable ll={}; +% t.traversal(&arrapp;list=ll); +% foreach i (ll) { +% print(i); +% } +% % print the minimum and maximum values +% ()=printf("min: %S - max: %S\n",t.min(),t.max()); +% +%\example +% % setup a tree for strings +% variable arr=["A","X","D","B","C"]; +% variable t=AVLTree(&strcmp); +% variable i; +% foreach i (arr) { +% t.insert(i); +% } +% t.print(); +% +% +%\seealso{AVLTree,sort} +%!%- +{ + variable T=@BinarySearchTree_Struct; + + if (_NARGS==1) { + variable cmpf=(); + T.key_compare=cmpf; + } + return T; +} + + +% +% AVL node specific functions and structures +% +private define new_AVLNode() { + variable T; + variable data=NULL; + if (_NARGS==1) { + T=(); + } + if (_NARGS==2) { + (T,data)=(); + } + + variable N=new_BinarySearchTreeNode(T,data); + return struct_combine(N,struct{height=0}); +} + +private define AVLNodeHeight(N) { + % get the node height + if (N==NULL) { + return -1; + } + return N.height; +} + +private define AVLNodeUpdateHeight(N) { + variable lh=AVLNodeHeight(N.child[TREE_LEFT]); + variable rh=AVLNodeHeight(N.child[TREE_RIGHT]); + if (lh>rh) { + N.height=lh+1; + } else { + N.height=rh+1; + } +} + +private define AVLNodeBalanceFactor(N) { + return AVLNodeHeight(N.child[TREE_RIGHT])-AVLNodeHeight(N.child[TREE_LEFT]); +} + +private define AVLTreeRotate(N,dir) { + % rotate tree below N by dir (where dir==TREE_LEFT or TREE_RIGHT) + variable sav=N.child[1-dir]; + + N.child[1-dir]=sav.child[dir]; + sav.child[dir]=N; + + AVLNodeUpdateHeight(N); + AVLNodeUpdateHeight(sav); + + return sav; +} + +private define AVLTreeRebalance(N) { + variable bf=AVLNodeBalanceFactor(N); + + % Left heavy? + if (bf < -1 ) { + if (AVLNodeBalanceFactor(N.child[TREE_LEFT])<=0) { + N=AVLTreeRotate(N,TREE_RIGHT); + } else { + N.child[TREE_LEFT]=AVLTreeRotate(N.child[TREE_LEFT],TREE_LEFT); + N=AVLTreeRotate(N,TREE_RIGHT); + } + + % if (AVLNodeBalanceFactor(N.child[TREE_LEFT])>0) { + % N.child[TREE_LEFT]=AVLTreeRotate(N.child[TREE_LEFT],TREE_LEFT); + % } + % N=AVLTreeRotate(N,TREE_RIGHT); + } + + % right heavy? + if (bf>1) { + if (AVLNodeBalanceFactor(N.child[TREE_RIGHT])>=0) { + N=AVLTreeRotate(N,TREE_LEFT); + } else { + N.child[TREE_RIGHT]=AVLTreeRotate(N.child[TREE_RIGHT],TREE_RIGHT); + N=AVLTreeRotate(N,TREE_LEFT); + } + % if (AVLNodeBalanceFactor(N.child[TREE_RIGHT])<0){ + % N.child[TREE_RIGHT]=AVLTreeRotate(N.child[TREE_RIGHT],TREE_RIGHT); + % } + % N=AVLTreeRotate(N,TREE_LEFT); + } + + return N; + +} + +private define AVLTreeInsert_helper(T,key,no); + +private define AVLTreeInsert_helper(T,key,no) { + % insert node N into binary tree no in tree T + if (no==NULL) { + return T.new_node(key); + } + + variable cmp=T.key_compare(key,no.data); + if (cmp<0) { + % key < no.data + no.child[TREE_LEFT]=AVLTreeInsert_helper(T,key,no.child[TREE_LEFT]); + } else { + if (cmp>0) { + no.child[TREE_RIGHT]=AVLTreeInsert_helper(T,key,no.child[TREE_RIGHT]); + } else { + % equal - not permitted, so throw an error + throw UsageError,sprintf("insert: Key %S already exists!\n",key); + } + } + + AVLNodeUpdateHeight(no); + + return AVLTreeRebalance(no); +} + +private define AVLTreeInsert (T,key) { + % insert a node containing key into the binary tree T + if (T.root==NULL) { + T.root=T.new_node(key); + return; + } + T.root=AVLTreeInsert_helper(T,key,T.root); +} + + +private define AVLTreeDeleteKey_helper(T,key,root) { + + variable N=BinarySearchTreeDeleteKey_helper_recursive(T,key,root); + if (N==NULL) { + return root; + } + AVLNodeUpdateHeight(N); + AVLTreeRebalance(N); +} + +private define AVLTreeDeleteKey(T,key) { + T.root=AVLTreeDeleteKey_helper(T,key,T.root); +} + +private define AVLTreeDeleteNode(T,N) { + T.root=AVLTreeDeleteKey_helper(T,N.key,T.root); +} + +private variable AVLTree_Struct=struct { + root=NULL, + insert=&AVLTreeInsert, + print=&BinarySearchTreePrint, + new_node=&new_AVLNode, + key_compare=&tree_key_compare, + search=&BinarySearchTreeSearch, + delete=&AVLTreeDeleteKey, + delete_node=&AVLTreeDeleteNode, + exists=&BinarySearchTreeKeyExists, + min=&BinarySearchTreeMinimum, + max=&BinarySearchTreeMaximum, + successor=&BinarySearchTreeNodeSuccessor, + predecessor=&BinarySearchTreeNodePredecessor, + traversal=&BinarySearchTraversal +}; + +define AVLTree() +%!%+ +%\function{AVLTree} +%\synopsis{Creates a new balanced binary search tree} +%\usage{Struct_Type=AVLTree();} +%\altusage{Struct_Type=AVLTree(cmpfunction);} +%\description +%A binary search tree is an object that permits fast +%searches and ordered operations on ordinal data types (i.e., +%data that can be sorted). +% +%An AVLTree is a balanced binary search tree, that is, a search +%tree with a structure that ensures that searches will always +%be close to O(log(N)), where N is the number of elements in the +%tree. +% +%The interface of AVLTree is identical to that of BinarySearchTree, +%see there for a description of the member functions. +% +%\seealso{BinarySearchTree} +%!%- +{ + variable T=@AVLTree_Struct; + + if (_NARGS==1) { + variable cmpf=(); + T.key_compare=cmpf; + } + return T; +} + + +% % +% % Red Black Tree specific functions +% % + +% private define RBRotateDirRoot(T,P,dir) { +% % T : red-black tree +% % P : root of subtree (may be the root of T) +% % dir: { TREE_LEFT, TREE_RIGHT } +% variable G=P.parent; +% variable S=P.child[1-dir]; + +% if (typeof(S)!=Struct_Type) { +% throw UsageError,"Need a true node!"; +% } + +% variable C = S.child[dir]; +% P.child[1-dir]=C; + +% if (typeof(C)==Struct_Type) { +% C.parent=P; +% } + +% S.child[dir]=P; +% P.parent=S; +% S.parent=G; + +% if (typeof(G)==Struct_Type) { +% G.child[ (P == G.child[TREE_RIGHT]) ? TREE_RIGHT : TREE_LEFT ] = S; +% } else { +% T.root=S; +% } +% return S; % new root of subtree +% } + +% define RBinsert1(T,N,P,dir) { +% % T: red-black tree +% % N: node to be inserted +% % P: parent node of N (may be NULL) +% % dir: side ( TREE_LEFT or TREE_RIGHT ) of P where to insert N + +% variable G; % parent node of P +% variable U; % uncle of N + +% N.color = RED; +% N.child[TREE_LEFT] = NULL; +% N.child[TREE_RIGHT] = NULL; +% N.parent = P; + +% if (typeof(P) != Struct_Type ) { % There is no parent +% T.root = N; % N is the new root of the tree T. +% return; % insertion complete +% } + +% P.child[dir]=N;% insert N as dir-child of P +% % start of the (do while)-loop: +% do { +% if (P.color == BLACK) { +% % Case_I1 (P black): +% return; % insertion complete +% } +% % From now on P is red. +% G=P.parent; +% if (typeof(G) != Struct_Type) { +% % P red and root +% P.color = BLACK; +% return; % insertion complete +% } +% % else: P red and G!=NULL. +% dir = childDir(P); % the side of parent G on which node P is located +% U = G.child[1-dir]; % uncle +% if (typeof(U) != Struct_Type || U.color == BLACK) { % considered black +% if (N == P.child[1-dir]) { +% % Case_I5 (P red && U black && N inner grandchild of G): +% RBRotateDirRoot(T,P,dir); % P is never the root +% N = P; % new current node +% P = G.child[dir]; % new parent of N +% % fall through to Case_I6 +% } +% %Case_I6 (P red && U black && N outer grandchild of G): +% RBRotateDirRoot(T,G,1-dir); % G may be the root +% P.color = BLACK; +% G.color = RED; +% return; % insertion complete +% } + +% % Case_I2 (P+U red): +% P.color = BLACK; +% U.color = BLACK; +% G.color = RED; +% N = G; % new current node +% % iterate 1 black level higher +% % (= 2 tree levels) +% P=N.parent; +% } while (typeof(P) == Struct_Type); +% % end of the (do while)-loop +% % Leaving the (do while)-loop (after having fallen through from Case_I2). + +% % Case_I3: N is the root and red. +% return; % insertion complete + +% } % end of RBinsert1 + +% define RBdelete(T,N) { +% % +% % delete node N from RB tree T +% % + +% % +% % simple root? +% % +% if (N==T.root) { +% if (typeof(N.child[TREE_LEFT])!=Struct_Type && typeof(N.child[TREE_RIGHT])!=Struct_Type) { +% % tree consists of one node only +% T.root=NULL; +% return; +% } + +% % now at least one child exists +% %%% if (typeof(XXX + +% } + +% } + +% define RBdelete2(T,N) { +% % +% % treatment of complex delete cases +% % +% % T red-black tree +% % N node to be deleted + +% variable P = N.parent; % -> parent node of N +% variable S; % sibling of N +% variable C; % close nephew +% variable D; % distant nephew + +% % P != NULL, since N is not the root. +% variable dir = childDir(N); % side of parent P on which the node N is located +% % Replace N at its parent P by NULL +% P.child[dir] = NULL; + +% % start of the (do while)-loop: +% do { +% dir = childDir(N); % side of parent P on which node N is located +% S = P.child[1-dir]; % sibling of N (has black height >= 1) +% D = S.child[1-dir]; % distant nephew +% C = S.child[ dir]; % close nephew +% if (S.color == RED) { +% % Case_D3: S red && P+C+D black: +% RBRotateDirRoot(T,P,dir); % P may be the root +% P.color = RED; +% S.color = BLACK; +% S = C; % != NIL +% % now: P red && S black +% D = S.child[1-dir]; % distant nephew +% if (typeof(D) == Struct_Type && D.color == RED) { +% % Case_D6: % D red && S black: +% RBRotateDirRoot(T,P,dir); % P may be the root +% S.color = P.color; +% P.color = BLACK; +% D.color = BLACK; +% return; +% } +% C = S.child[ dir]; % close nephew +% if (typeof(C) == Struct_Type && C.color == RED) { +% % Case_D5: % C red && S+D black: +% RBRotateDirRoot(T,S,1-dir); % S is never the root +% S.color = RED; +% C.color = BLACK; +% D = S; +% S = C; +% RBRotateDirRoot(T,P,dir); % P may be the root +% S.color = P.color; +% P.color = BLACK; +% D.color = BLACK; +% return; +% } +% % Otherwise C+D considered black. +% S.color = RED; +% P.color = BLACK; +% return; +% } +% % S is black: +% if (typeof(D) == Struct_Type && D.color == RED) { % not considered black +% % Case_D6: % D red && S black: +% RBRotateDirRoot(T,P,dir); % P may be the root +% S.color = P.color; +% P.color = BLACK; +% D.color = BLACK; +% return; % deletion complete +% } +% if (typeof(C) == Struct_Type && C.color == RED) {% not considered black +% % C red && S+D black: +% RBRotateDirRoot(T,S,1-dir); % S is never the root +% S.color = RED; +% C.color = BLACK; +% D = S; +% S = C; +% RBRotateDirRoot(T,P,dir); % P may be the root +% S.color = P.color; +% P.color = BLACK; +% D.color = BLACK; +% return; +% } +% % Here both nephews are == NULL (first iteration) or black (later). +% if (P.color == RED) { +% S.color = RED; +% P.color = BLACK; +% return; +% } +% % Case_D1 (P+C+S+D black): +% S.color = RED; +% N = P; % new current node (maybe the root) +% % iterate 1 black level +% % (= 1 tree level) higher +% P=N.parent; +% } while (typeof(P) == Struct_Type); +% % end of the (do while)-loop +% return; + +% } % end of RBdelete + +% private define RedBlackTreeInsert(T,key) { +% % insert a node containing key into RB tree T + +% variable N=new_RBNode(data); + +% if (typeof(T.root)!=Struct_Type) { +% T.root=N; +% T.root.color=BLACK; +% return; +% } + +% variable x=T.root; + +% % find the future parent of our node +% while (typeof(x) == Struct_Type ) { +% variable y=x; % +% if (N.data < x.data) { +% x=x.child[TREE_LEFT]; +% } else { +% x=x.child[TREE_RIGHT]; +% } +% } +% N.parent=y; + +% % where to insert +% variable dir=(T.key_compare(N.data,y.data)==-1) ? TREE_LEFT : TREE_RIGHT; + +% RBinsert1(T,N,N.parent,dir); +% } + + + +% variable RedBlackTree_Struct=struct { +% root=NULL, +% new_node=&new_RedBlackTreeNode, +% insert=&RedBlackTreeInsert, +% print=&BinarySearchTreePrint, +% key_compare=&tree_key_compare, +% search=&BinarySearchTreeSearch, +% delete=&RBdelete, +% exists=&BinarySearchTreeExists, +% min=&BinarySearchTreeMinimum, +% max=&BinarySearchTreeMaximum, +% successor=&BinaryNodeSuccessor, +% predecessor=&BinaryNodePredecessor +% }; + + + +% private define new_RBTree() { +% return @RedBlackTree_Struct; +% } + + + +% % function joinRightRB(TL, k, TR): +% % if (TL.color=black) and (TL.blackHeight=TR.blackHeight): +% % return Node(TL, {k,red},TR) +% % T'=Node(TL.left, {TL.key,TL.color},joinRightRB(TL.right,k,TR)) +% % if (TL.color=black) and (T'.right.color=T'.right.right.color=red): +% % T'.right.right.color=black; +% % return rotateLeft(T',N,T') +% % return T' /* T''[recte T'] */ + +% % function joinLeftRB(TL, k, TR): +% % /* symmetric to joinRightRB */ + +% % function join(TL, k, TR): +% % if TL.blackHeight>TR.blackHeight: +% % T'=joinRightRB(TL,k,TR) +% % if (T'.color=red) and (T'.right.color=red): +% % T'.color=black +% % return T' +% % if TR.blackHeight>TL.blackHeight: +% % /* symmetric */ +% % if (TL.color=black) and (TR.color=black): +% % return Node(TL,{k,red},TR) +% % return Node(TL,{k,blac},TR) + +% % function split(T, k): +% % if (T = nil) return (nil, false, nil) +% % if (k = T.key) return (T.left, true, T.right) +% % if (k < T.key): +% % (L',b,R') = split(T.left, k) +% % return (L',b,join(R',T.key,T.right)) +% % (L',b,R') = split(T.right, k) +% % return (join(T.left,T.key,L'),b,T.right) + +% % function union(t1, t2): +% % if t1 = nil return t2 +% % if t2 = nil return t1 +% % (L1,b,R1)=split(t1,t2.key) +% % proc1=start: +% % TL=union(L1,t2.left) +% % proc2=start: +% % TR=union(R1,t2.right) +% % wait all proc1,proc2 +% % return join(TL, t2.key, TR) -- GitLab From 285119a499d46eb7500054ccff8820b43656cc9d Mon Sep 17 00:00:00 2001 From: Joern Wilms Date: Wed, 24 Aug 2022 22:18:50 +0200 Subject: [PATCH 7/9] added clipping of polylines with a polygon --- src/slang/geometry/clip_polyline_polygon.sl | 183 ++++++++++++++++++++ 1 file changed, 183 insertions(+) create mode 100755 src/slang/geometry/clip_polyline_polygon.sl diff --git a/src/slang/geometry/clip_polyline_polygon.sl b/src/slang/geometry/clip_polyline_polygon.sl new file mode 100755 index 00000000..e4589edd --- /dev/null +++ b/src/slang/geometry/clip_polyline_polygon.sl @@ -0,0 +1,183 @@ +define clip_polyline_polygon(polyline,polygon) +%!%+ +%\function{clip_polyline_polygon} +%\usage{clipped=clip_polyline_polygon(line,polygon) +%\qualifiers{ +%\qualifier{evenodd}{: use the crossing number method} +%\qualifier{crossing}{: use the crossing number method} +%\qualifier{winding}{: use the winding number method (the default)} +%} +%\description +%This function clips a polyline defined by line=struct{x=[],y=[]}, +%where x and y are the points connected by the polyline, against +%a closed polygon=struct{x=[],y=[]}, where again x and y are the points +%connected by the polygon and where x[0]==x[-1] and y[0]==y[-1]. +% +%The function returns a list of polylines that contain the segments +%of the line that are inside of the polygon. This list can be empty if +%there is no overlap. +% +%For complex polygons with intersecting segments, the qualifier defines +%what constitues the inside of the polygon. +% +%If you want to clip against a simple rectangle, use clip_polyline_rectangle +%for a faster algorithm. +% +%\seealso{point_in_polygon,clip_polyline_rectangle,clip_points_polygon} +%!%- +{ + variable evenodd=(qualifier_exists("crossing") || qualifier_exists("evenodd")); + + + % + % intersect each segment of poly with all segments of + % the polygon + % FIXME: this is ok for polygons and polylines with a small number + % of segments. For polygons with a large number of segments we still + % need to implement a version of the function that is based on + % a binary search tree. + % + variable i,j; + + variable outpolyx={}; + variable outpolyy={}; + variable state={}; + + variable crossing=2; + variable oldpt=1; + + variable npolygon=length(polygon.x); + variable nline=length(polyline.x); + % + % search for intersections of the segments of the polyline with + % all segments of the polygon (so this can be expensive) + % + % we generate a new temporary polyline, outpoly, that contains + % the points of the polyline and all crossing points, + % remembering whether we're dealing with a new crossing point or + % not + % + % FIXME: this can be memory intensive, since we're temporarily + % duplicating polyline. A better approach might be to remember + % the alphalists for each polyline segment (see below) and then + % use these in the polyline splitting step that follows. This is a + % bit more complicated to implement and up to now I have not + % encountered real memory problems that have forced me to takle + % this + % + _for i(0,nline-2) { + list_append(outpolyx,polyline.x[i]); + list_append(outpolyy,polyline.y[i]); + list_append(state,oldpt); + + variable alphalist={}; + _for j(0,npolygon-2) { + variable alphap,alphaq; + variable res=line_intersect(polyline.x[i],polyline.y[i],polyline.x[i+1],polyline.y[i+1], + polygon.x[j],polygon.y[j],polygon.x[j+1],polygon.y[j+1],&alphap,&alphaq); + if (res==1) { + list_append(alphalist,alphap); + } + } + + % if there are intersections: add them in ordered form (i.e., in the + % direction of the polyline) to the outpoly + if (length(alphalist)>0) { + alphalist=list_to_array(alphalist); + variable ndx=array_sort(alphalist); + foreach alphap ( alphalist[ndx] ) { + variable xx=polyline.x[i]+alphap*(polyline.x[i+1]-polyline.x[i]); + variable yy=polyline.y[i]+alphap*(polyline.y[i+1]-polyline.y[i]); + list_append(outpolyx,xx); + list_append(outpolyy,yy); + list_append(state,crossing); + } + } + } + + % don't forget the last point + list_append(outpolyx,polyline.x[-1]); + list_append(outpolyy,polyline.y[-1]); + list_append(state,oldpt); + + outpolyx=list_to_array(outpolyx); + outpolyy=list_to_array(outpolyy); + + % + % second step: walk along the polyline and cut it apart. + % + % This is easy in the case of the even odd rule, but a bit + % more complex in the case of the winding number/crossing method + % + + % determine whether the first point is inside or outside + % (__qualifiers given to ensure that evenodd is taken into account) + % this gives the "state" of the current segment + variable inside=point_in_polygon(outpolyx[0],outpolyy[0],polygon.x,polygon.y;;__qualifiers()); + + % this is where the individual segments will end up + variable polypieces={}; + + % coordinates of the current segment + variable segx={}; + variable segy={}; + + _for i(0,length(outpolyx)-2) { + if (inside) { + list_append(segx,outpolyx[i]); + list_append(segy,outpolyy[i]); + } + + variable newseg=0; % 1 if we start a new segment + if (state[i]==crossing) { + if (evenodd) { + % even odd rule: the state of the segment changes with each crossing + newseg=1; + } else { + % winding number rule: it is possible to either remain inside or go outside. + % Test this by looking at the location of the middle point between point i and i+1, + % to avoid round (the point in the middle is guaranteed not to be on a crossing, + % after all) + % + variable xm=(outpolyx[i]+outpolyx[i+1])/2.; + variable ym=(outpolyy[i]+outpolyy[i+1])/2.; + % no need for qualifiers: the default of point_in_polygin is the + % winding number/crossing number + variable segpt=point_in_polygon(xm,ym,polygon.x,polygon.y); + % a new segment starts if the state changes + % (remember that inside=0 is outside) + newseg=(segpt!=inside); + } + + % + % if we're starting a new segment, save the polygon or start a new one + % + if (newseg) { + if (inside) { + % we're leaving the polygon, so save the current + % segment and start a new one + list_append(polypieces,struct{x=list_to_array(segx),y=list_to_array(segy)}); + inside=0; + segx={}; + segy={}; + } else { + % we're entering the polygon + inside=1; + list_append(segx,outpolyx[i]); + list_append(segy,outpolyy[i]); + } + } + } + } + + % the last point is never a crossing point, so the state can't change + if (inside) { + list_append(segx,outpolyx[-1]); + list_append(segy,outpolyy[-1]); + list_append(polypieces,struct{x=list_to_array(segx),y=list_to_array(segy)}); + } + + % and we're done + return polypieces; +} + -- GitLab From 06861ac8eca6e2d9d140962903bdf24f7173c580 Mon Sep 17 00:00:00 2001 From: Joern Wilms Date: Wed, 24 Aug 2022 22:22:57 +0200 Subject: [PATCH 8/9] prevented namespace pollution --- src/slang/trees/avltree.sl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/slang/trees/avltree.sl b/src/slang/trees/avltree.sl index 8f4af529..e6dbcc39 100644 --- a/src/slang/trees/avltree.sl +++ b/src/slang/trees/avltree.sl @@ -106,7 +106,7 @@ private define BinarySearchTreeMinimum(N) { return N.data; } -define BinarySearchTreeMaximum(N) { +private define BinarySearchTreeMaximum(N) { % return the largest key in a binary search tree % if qualifier "node" is set returns the node, otherwise the key @@ -125,8 +125,8 @@ define BinarySearchTreeMaximum(N) { } % iterative deletion of a key in a binary search tree -define BinarySearchTreeDeleteKey_helper_iterative(T,key,N); -define BinarySearchTreeDeleteKey_helper_iterative(T,key,N) { +private define BinarySearchTreeDeleteKey_helper_iterative(T,key,N); +private define BinarySearchTreeDeleteKey_helper_iterative(T,key,N) { % % deletion of node with key K from binary tree (version 1) % @@ -187,9 +187,9 @@ define BinarySearchTreeDeleteKey_helper_iterative(T,key,N) { % % full recursive tree deletion (needed for the AVL tree) % -define BinarySearchTreeDeleteKey_helper_recursive(T,key,N); +private define BinarySearchTreeDeleteKey_helper_recursive(T,key,N); -define BinarySearchTreeDeleteKey_helper_recursive(T,key,N) { +private define BinarySearchTreeDeleteKey_helper_recursive(T,key,N) { % no node at current position, end recursion if (N==NULL) { return NULL; @@ -234,17 +234,17 @@ define BinarySearchTreeDeleteKey_helper_recursive(T,key,N) { return N; } -define BinarySearchTreeDeleteKey(T,key) { +private define BinarySearchTreeDeleteKey(T,key) { ()=BinarySearchTreeDeleteKey_helper_iterative(T,key,T.root); } -define BinarySearchTreeDeleteNode(T,N) { +private define BinarySearchTreeDeleteNode(T,N) { ()=BinarySearchTreeDeleteKey_helper_iterative(T,N.key,T.root); } % find node for which data==key % returns the Node of Null, if the key does not exist -define BinarySearchTreeSearch(T,key) { +private define BinarySearchTreeSearch(T,key) { variable N=T.root; while (N!=NULL && N.data != key ) { if (key < N.data) { @@ -368,7 +368,7 @@ private define BinarySearchTreePrint(T) { BinarySearchTreeNodePrint(T.root,""); } -variable BinarySearchTree_Struct=struct { +private variable BinarySearchTree_Struct=struct { root=NULL, insert=&BinarySearchTreeInsert, print=&BinarySearchTreePrint, -- GitLab From 4764616801a05354bb4a4fad286187c51a8825dd Mon Sep 17 00:00:00 2001 From: Joern Wilms Date: Fri, 26 Aug 2022 01:03:30 +0200 Subject: [PATCH 9/9] added fast function to ceck for overlapping rectangles --- src/slang/geometry/rectangles_overlap.sl | 46 ++++++++++++++++++++++++ 1 file changed, 46 insertions(+) create mode 100644 src/slang/geometry/rectangles_overlap.sl diff --git a/src/slang/geometry/rectangles_overlap.sl b/src/slang/geometry/rectangles_overlap.sl new file mode 100644 index 00000000..0b0cd5ec --- /dev/null +++ b/src/slang/geometry/rectangles_overlap.sl @@ -0,0 +1,46 @@ +define rectangles_overlap (src,clp) +%!%+ +%\function{rectangles_overlap} +%\synopsis{Checks whether two rectangles overlap} +%\usage{res=rectangles_overlap(src,clp)} +%\description +% +%This function checks whether the two rectangles src and clp +%overlap (return value=1) or not (return value =0). +% +% The rectangles are defined by (xmin,ymin,xmax,ymax) +% either as a 4 element array in this order or as a +% struct with these tags. Both boxes must be defined +% using the same format. +% +%!%- +{ + if (typeof(src)!=typeof(clp)) { + throw UsageError,"%s: src and clp must have the same type!",_function_name(); + } + + variable xmin,xmax,ymin,ymax; + if (typeof(src)==Array_Type) { + % (sxmax<=cxmin || sxmin>=cxmax) + if (src[2]<=clp[0] || src[0]>=clp[2]) { + return 0; + } + + % if (symax<=cymin || symin>=cymax) + if (src[3]<=clp[1] || src[1]>=clp[3]) { + return 0; + } + return 1; + } + + if (src.xmax<=clp.xmin || src.xmin>=clp.xmax) { + return 0; + } + + if (src.ymax<=clp.ymin || src.ymin>=clp.ymax) { + return 0; + } + + return 1; + +} -- GitLab