From: John Darrington Date: Thu, 23 Jul 2009 05:05:41 +0000 (+0200) Subject: Merge commit 'origin/roc' X-Git-Tag: build37~50 X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?p=pspp-builds.git;a=commitdiff_plain;h=b40baf410822471fbdeeec553693619d60d7c7b6;hp=231ebe7e6d5c8c8eba47560793b5fa53c013b87f Merge commit 'origin/roc' --- diff --git a/doc/statistics.texi b/doc/statistics.texi index 921ea857..76925d7e 100644 --- a/doc/statistics.texi +++ b/doc/statistics.texi @@ -15,6 +15,7 @@ far. * RANK:: Compute rank scores. * REGRESSION:: Linear regression. * RELIABILITY:: Reliability analysis. +* ROC:: Receiver Operating Characteristic. @end menu @node DESCRIPTIVES @@ -950,3 +951,73 @@ analysis tested against the totals. +@node ROC +@section ROC + +@vindex ROC +@cindex Receiver Operating Characterstic +@cindex Area under curve + +@display +ROC @var{var_list} BY @var{state_var} (@var{state_value}) + /PLOT = @{ CURVE [(REFERENCE)], NONE @} + /PRINT = [ SE ] [ COORDINATES ] + /CRITERIA = [ CUTOFF(@{INCLUDE,EXCLUDE@}) ] + [ TESTPOS (@{LARGE,SMALL@}) ] + [ CI (@var{confidence}) ] + [ DISTRIBUTION (@{FREE, NEGEXPO @}) ] + /MISSING=@{EXCLUDE,INCLUDE@} +@end display + + +The @cmd{ROC} command is used to plot the receiver operating characteristic curve +of a dataset, and to estimate the area under the curve. +This is useful for analysing the efficacy of a variable as a predictor of a state of nature. + +The mandatory @var{var_list} is the list of predictor variables. +The variable @var{state_var} is the variable whose values represent the actual states, +and @var{state_value} is the value of this variable which represents the positive state. + +The optional subcommand PLOT is used to determine if and how the ROC curve is drawn. +The keyword CURVE means that the ROC curve should be drawn, and the optional keyword REFERENCE, +which should be enclosed in parentheses, says that the diagonal reference line should be drawn. +If the keyword NONE is given, then no ROC curve is drawn. +By default, the curve is drawn with no reference line. + +The optional subcommand PRINT determines which additional tables should be printed. +Two additional tables are available. +The SE keyword says that standard error of the area under the curve should be printed as well as +the area itself. +In addition, a p-value under the null hypothesis that the area under the curve equals 0.5 will be +printed. +The COORDINATES keyword says that a table of coordinates of the ROC curve should be printed. + +The CRITERIA subcommand has four optional parameters: +@itemize @bullet +@item The TESTPOS parameter may be LARGE or SMALL. +LARGE is the default, and says that larger values in the predictor variables are to be +considered positive. SMALL indicates that smaller values should be considered positive. + +@item The CI parameter specifies the confidence interval that should be printed. +It has no effect if the SE keyword in the PRINT subcommand has not been given. + +@item The DISTRIBUTION parameter determines the method to be used when estimating the area +under the curve. +There are two possibilities, @i{viz}: FREE and NEGEXPO. +The FREE method uses a non-parametric estimate, and the NEGEXPO method a bi-negative +exponential distribution estimate. +The NEGEXPO method should only be used when the number of positive actual states is +equal to the number of negative actual states. +The default is FREE. + +@item The CUTOFF parameter is for compatibility and is ignored. +@end itemize + +The MISSING subcommand determines whether user missing values are to +be included or excluded in the analysis. The default behaviour is to +exclude them. +Cases are excluded on a listwise basis; if any of the variables in @var{var_list} +or if the variable @var{state_var} is missing, then the entire case will be +excluded. + + diff --git a/src/data/casereader-translator.c b/src/data/casereader-translator.c index 548c22fb..feffa15c 100644 --- a/src/data/casereader-translator.c +++ b/src/data/casereader-translator.c @@ -19,6 +19,7 @@ #include #include +#include #include #include @@ -359,3 +360,132 @@ car_translate (struct ccase *input, void *car_) } + + +struct consolidator +{ + const struct variable *key; + const struct variable *weight; + double cc; + double prev_cc; + + casenumber n; + struct casereader *clone; + struct caseproto *proto; + int direction; +}; + +static bool +uniquify (const struct ccase *c, void *aux) +{ + struct consolidator *cdr = aux; + const union value *current_value = case_data (c, cdr->key); + const int key_width = var_get_width (cdr->key); + const double weight = cdr->weight ? case_data (c, cdr->weight)->f : 1.0; + const struct ccase *next_case = casereader_peek (cdr->clone, cdr->n + 1); + int dir = 0; + + cdr->n ++; + cdr->cc += weight; + + if ( NULL == next_case) + goto end; + + dir = value_compare_3way (case_data (next_case, cdr->key), + current_value, key_width); + if ( dir != 0 ) + { + /* Insist that the data are sorted */ + assert (cdr->direction == 0 || dir == cdr->direction); + cdr->direction = dir; + goto end; + } + + return false; + + end: + cdr->prev_cc = cdr->cc; + cdr->cc = 0; + return true; +} + + + +static struct ccase * +consolodate_weight (struct ccase *input, void *aux) +{ + struct consolidator *cdr = aux; + struct ccase *c; + + if (cdr->weight) + { + c = case_unshare (input); + case_data_rw (c, cdr->weight)->f = cdr->prev_cc; + } + else + { + c = case_unshare_and_resize (input, cdr->proto); + case_data_rw_idx (c, caseproto_get_n_widths (cdr->proto) - 1)->f = cdr->prev_cc; + } + + return c; +} + + +static bool +uniquify_destroy (void *aux) +{ + struct consolidator *cdr = aux; + + casereader_destroy (cdr->clone); + caseproto_unref (cdr->proto); + free (cdr); + + return true; +} + + + +/* Returns a new casereader which is based upon INPUT, but which contains a maximum + of one case for each distinct value of KEY. + If WEIGHT is non-null, then the new casereader's values for this variable + will be the sum of all values matching KEY. + IF WEIGHT is null, then the new casereader will have an additional numeric + value appended, which will contain the total number of cases containing + KEY. + INPUT must be sorted on KEY +*/ +struct casereader * +casereader_create_distinct (struct casereader *input, + const struct variable *key, + const struct variable *weight) +{ + struct casereader *u ; + struct casereader *ud ; + struct caseproto *output_proto = caseproto_ref (casereader_get_proto (input)); + + struct consolidator *cdr = xmalloc (sizeof (*cdr)); + cdr->n = 0; + cdr->key = key; + cdr->weight = weight; + cdr->cc = 0; + cdr->clone = casereader_clone (input); + cdr->direction = 0; + + if ( NULL == cdr->weight ) + output_proto = caseproto_add_width (output_proto, 0); + + cdr->proto = output_proto; + + u = casereader_create_filter_func (input, uniquify, + NULL, cdr, NULL); + + ud = casereader_create_translator (u, + output_proto, + consolodate_weight, + uniquify_destroy, + cdr); + + return ud; +} + diff --git a/src/data/casereader.h b/src/data/casereader.h index d4f2966a..3b903bbf 100644 --- a/src/data/casereader.h +++ b/src/data/casereader.h @@ -140,5 +140,10 @@ casereader_create_append_rank (struct casereader *, enum rank_error *err, distinct_func *distinct_callback, void *aux); +struct casereader * +casereader_create_distinct (struct casereader *input, + const struct variable *key, + const struct variable *weight); + #endif /* data/casereader.h */ diff --git a/src/data/subcase.c b/src/data/subcase.c index be586096..6ffaa4c2 100644 --- a/src/data/subcase.c +++ b/src/data/subcase.c @@ -64,6 +64,16 @@ subcase_init_var (struct subcase *sc, const struct variable *var, subcase_add_var (sc, var, direction); } + +void +subcase_init (struct subcase *sc, int index, int width, + enum subcase_direction direction) +{ + subcase_init_empty (sc); + subcase_add (sc, index, width, direction); +} + + /* Removes all the fields from SC. */ void subcase_clear (struct subcase *sc) @@ -89,6 +99,7 @@ subcase_destroy (struct subcase *sc) caseproto_unref (sc->proto); } + /* Add a field for VAR to SC, with DIRECTION as the sort order. Returns true if successful, false if VAR already has a field in SC. */ @@ -96,7 +107,17 @@ bool subcase_add_var (struct subcase *sc, const struct variable *var, enum subcase_direction direction) { - size_t case_index = var_get_case_index (var); + return subcase_add (sc, var_get_case_index (var), + var_get_width (var), direction); +} + +/* Add a field for CASE_INDEX, WIDTH to SC, with DIRECTION as the sort order. + Returns true if successful, false if CASE_INDEX already has a field + in SC. */ +bool +subcase_add (struct subcase *sc, int case_index, int width, + enum subcase_direction direction) +{ struct subcase_field *field; size_t i; @@ -107,7 +128,7 @@ subcase_add_var (struct subcase *sc, const struct variable *var, sc->fields = xnrealloc (sc->fields, sc->n_fields + 1, sizeof *sc->fields); field = &sc->fields[sc->n_fields++]; field->case_index = case_index; - field->width = var_get_width (var); + field->width = width; field->direction = direction; invalidate_proto (sc); return true; diff --git a/src/data/subcase.h b/src/data/subcase.h index 050cf17d..6e59da1d 100644 --- a/src/data/subcase.h +++ b/src/data/subcase.h @@ -53,10 +53,16 @@ void subcase_init_vars (struct subcase *, const struct variable *const *, size_t n_vars); void subcase_init_var (struct subcase *, const struct variable *, enum subcase_direction); +void subcase_init (struct subcase *, int index, int width, + enum subcase_direction); + void subcase_clone (struct subcase *, const struct subcase *); void subcase_clear (struct subcase *); void subcase_destroy (struct subcase *); +bool subcase_add (struct subcase *sc, int index, int width, + enum subcase_direction direction); + bool subcase_add_var (struct subcase *, const struct variable *, enum subcase_direction); diff --git a/src/language/command.def b/src/language/command.def index 4c8df335..fa1bb1e8 100644 --- a/src/language/command.def +++ b/src/language/command.def @@ -117,6 +117,7 @@ DEF_CMD (S_DATA, 0, "RANK", cmd_rank) DEF_CMD (S_DATA, 0, "REGRESSION", cmd_regression) DEF_CMD (S_DATA, 0, "RELIABILITY", cmd_reliability) DEF_CMD (S_DATA, 0, "RENAME VARIABLES", cmd_rename_variables) +DEF_CMD (S_DATA, 0, "ROC", cmd_roc) DEF_CMD (S_DATA, 0, "SAMPLE", cmd_sample) DEF_CMD (S_DATA, 0, "SAVE", cmd_save) DEF_CMD (S_DATA, 0, "SORT CASES", cmd_sort_cases) @@ -237,7 +238,6 @@ UNIMPL_CMD ("REPEATING DATA", "Specify multiple cases per input record") UNIMPL_CMD ("REPORT", "Pretty print working file") UNIMPL_CMD ("RESTORE", "Restore settings") UNIMPL_CMD ("RMV", "Replace missing values") -UNIMPL_CMD ("ROC", "Receiver operating characteristic") UNIMPL_CMD ("SAVE TRANSLATE", "Save to foriegn format") UNIMPL_CMD ("SCRIPT", "Run script file") UNIMPL_CMD ("SEASON", "Estimate seasonal factors") diff --git a/src/language/stats/automake.mk b/src/language/stats/automake.mk index 5aee445c..1a68b906 100644 --- a/src/language/stats/automake.mk +++ b/src/language/stats/automake.mk @@ -32,10 +32,12 @@ language_stats_sources = \ src/language/stats/freq.h \ src/language/stats/npar-summary.c \ src/language/stats/npar-summary.h \ - src/language/stats/wilcoxon.c \ - src/language/stats/wilcoxon.h \ + src/language/stats/roc.c \ + src/language/stats/roc.h \ src/language/stats/sign.c \ - src/language/stats/sign.h + src/language/stats/sign.h \ + src/language/stats/wilcoxon.c \ + src/language/stats/wilcoxon.h all_q_sources += $(src_language_stats_built_sources:.c=.q) EXTRA_DIST += $(src_language_stats_built_sources:.c=.q) diff --git a/src/language/stats/roc.c b/src/language/stats/roc.c new file mode 100644 index 00000000..4cd593a2 --- /dev/null +++ b/src/language/stats/roc.c @@ -0,0 +1,1227 @@ +/* PSPP - a program for statistical analysis. + Copyright (C) 2009 Free Software Foundation, Inc. + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . */ + +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + + +#include + +#include +#include + +#include +#include + +#include "gettext.h" +#define _(msgid) gettext (msgid) +#define N_(msgid) msgid + +struct cmd_roc +{ + size_t n_vars; + const struct variable **vars; + const struct dictionary *dict; + + const struct variable *state_var ; + union value state_value; + + /* Plot the roc curve */ + bool curve; + /* Plot the reference line */ + bool reference; + + double ci; + + bool print_coords; + bool print_se; + bool bi_neg_exp; /* True iff the bi-negative exponential critieria + should be used */ + enum mv_class exclude; + + bool invert ; /* True iff a smaller test result variable indicates + a positive result */ + + double pos; + double neg; + double pos_weighted; + double neg_weighted; +}; + +static int run_roc (struct dataset *ds, struct cmd_roc *roc); + +int +cmd_roc (struct lexer *lexer, struct dataset *ds) +{ + struct cmd_roc roc ; + const struct dictionary *dict = dataset_dict (ds); + + roc.vars = NULL; + roc.n_vars = 0; + roc.print_se = false; + roc.print_coords = false; + roc.exclude = MV_ANY; + roc.curve = true; + roc.reference = false; + roc.ci = 95; + roc.bi_neg_exp = false; + roc.invert = false; + roc.pos = roc.pos_weighted = 0; + roc.neg = roc.neg_weighted = 0; + roc.dict = dataset_dict (ds); + + if (!parse_variables_const (lexer, dict, &roc.vars, &roc.n_vars, + PV_APPEND | PV_NO_DUPLICATE | PV_NUMERIC)) + goto error; + + if ( ! lex_force_match (lexer, T_BY)) + { + goto error; + } + + roc.state_var = parse_variable (lexer, dict); + + if ( !lex_force_match (lexer, '(')) + { + goto error; + } + + parse_value (lexer, &roc.state_value, var_get_width (roc.state_var)); + + + if ( !lex_force_match (lexer, ')')) + { + goto error; + } + + + while (lex_token (lexer) != '.') + { + lex_match (lexer, '/'); + if (lex_match_id (lexer, "MISSING")) + { + lex_match (lexer, '='); + while (lex_token (lexer) != '.' && lex_token (lexer) != '/') + { + if (lex_match_id (lexer, "INCLUDE")) + { + roc.exclude = MV_SYSTEM; + } + else if (lex_match_id (lexer, "EXCLUDE")) + { + roc.exclude = MV_ANY; + } + else + { + lex_error (lexer, NULL); + goto error; + } + } + } + else if (lex_match_id (lexer, "PLOT")) + { + lex_match (lexer, '='); + if (lex_match_id (lexer, "CURVE")) + { + roc.curve = true; + if (lex_match (lexer, '(')) + { + roc.reference = true; + lex_force_match_id (lexer, "REFERENCE"); + lex_force_match (lexer, ')'); + } + } + else if (lex_match_id (lexer, "NONE")) + { + roc.curve = false; + } + else + { + lex_error (lexer, NULL); + goto error; + } + } + else if (lex_match_id (lexer, "PRINT")) + { + lex_match (lexer, '='); + while (lex_token (lexer) != '.' && lex_token (lexer) != '/') + { + if (lex_match_id (lexer, "SE")) + { + roc.print_se = true; + } + else if (lex_match_id (lexer, "COORDINATES")) + { + roc.print_coords = true; + } + else + { + lex_error (lexer, NULL); + goto error; + } + } + } + else if (lex_match_id (lexer, "CRITERIA")) + { + lex_match (lexer, '='); + while (lex_token (lexer) != '.' && lex_token (lexer) != '/') + { + if (lex_match_id (lexer, "CUTOFF")) + { + lex_force_match (lexer, '('); + if (lex_match_id (lexer, "INCLUDE")) + { + roc.exclude = MV_SYSTEM; + } + else if (lex_match_id (lexer, "EXCLUDE")) + { + roc.exclude = MV_USER | MV_SYSTEM; + } + else + { + lex_error (lexer, NULL); + goto error; + } + lex_force_match (lexer, ')'); + } + else if (lex_match_id (lexer, "TESTPOS")) + { + lex_force_match (lexer, '('); + if (lex_match_id (lexer, "LARGE")) + { + roc.invert = false; + } + else if (lex_match_id (lexer, "SMALL")) + { + roc.invert = true; + } + else + { + lex_error (lexer, NULL); + goto error; + } + lex_force_match (lexer, ')'); + } + else if (lex_match_id (lexer, "CI")) + { + lex_force_match (lexer, '('); + lex_force_num (lexer); + roc.ci = lex_number (lexer); + lex_get (lexer); + lex_force_match (lexer, ')'); + } + else if (lex_match_id (lexer, "DISTRIBUTION")) + { + lex_force_match (lexer, '('); + if (lex_match_id (lexer, "FREE")) + { + roc.bi_neg_exp = false; + } + else if (lex_match_id (lexer, "NEGEXPO")) + { + roc.bi_neg_exp = true; + } + else + { + lex_error (lexer, NULL); + goto error; + } + lex_force_match (lexer, ')'); + } + else + { + lex_error (lexer, NULL); + goto error; + } + } + } + else + { + lex_error (lexer, NULL); + break; + } + } + + if ( ! run_roc (ds, &roc)) + goto error; + + free (roc.vars); + return CMD_SUCCESS; + + error: + free (roc.vars); + return CMD_FAILURE; +} + + + + +static void +do_roc (struct cmd_roc *roc, struct casereader *group, struct dictionary *dict); + + +static int +run_roc (struct dataset *ds, struct cmd_roc *roc) +{ + struct dictionary *dict = dataset_dict (ds); + bool ok; + struct casereader *group; + + struct casegrouper *grouper = casegrouper_create_splits (proc_open (ds), dict); + while (casegrouper_get_next_group (grouper, &group)) + { + do_roc (roc, group, dataset_dict (ds)); + } + ok = casegrouper_destroy (grouper); + ok = proc_commit (ds) && ok; + + return ok; +} + +#if 0 +static void +dump_casereader (struct casereader *reader) +{ + struct ccase *c; + struct casereader *r = casereader_clone (reader); + + for ( ; (c = casereader_read (r) ); case_unref (c)) + { + int i; + for (i = 0 ; i < case_get_value_cnt (c); ++i) + { + printf ("%g ", case_data_idx (c, i)->f); + } + printf ("\n"); + } + + casereader_destroy (r); +} +#endif + + +/* + Return true iff the state variable indicates that C has positive actual state. + + As a side effect, this function also accumulates the roc->{pos,neg} and + roc->{pos,neg}_weighted counts. + */ +static bool +match_positives (const struct ccase *c, void *aux) +{ + struct cmd_roc *roc = aux; + const struct variable *wv = dict_get_weight (roc->dict); + const double weight = wv ? case_data (c, wv)->f : 1.0; + + const bool positive = + ( 0 == value_compare_3way (case_data (c, roc->state_var), &roc->state_value, + var_get_width (roc->state_var))); + + if ( positive ) + { + roc->pos++; + roc->pos_weighted += weight; + } + else + { + roc->neg++; + roc->neg_weighted += weight; + } + + return positive; +} + + +#define VALUE 0 +#define N_EQ 1 +#define N_PRED 2 + +/* Some intermediate state for calculating the cutpoints and the + standard error values */ +struct roc_state +{ + double auc; /* Area under the curve */ + + double n1; /* total weight of positives */ + double n2; /* total weight of negatives */ + + /* intermediates for standard error */ + double q1hat; + double q2hat; + + /* intermediates for cutpoints */ + struct casewriter *cutpoint_wtr; + struct casereader *cutpoint_rdr; + double prev_result; + double min; + double max; +}; + +#define CUTPOINT 0 +#define TP 1 +#define FN 2 +#define TN 3 +#define FP 4 + + +/* + Return a new casereader based upon CUTPOINT_RDR. + The number of "positive" cases are placed into + the position TRUE_INDEX, and the number of "negative" cases + into FALSE_INDEX. + POS_COND and RESULT determine the semantics of what is + "positive". + WEIGHT is the value of a single count. + */ +static struct casereader * +accumulate_counts (struct casereader *cutpoint_rdr, + double result, double weight, + bool (*pos_cond) (double, double), + int true_index, int false_index) +{ + const struct caseproto *proto = casereader_get_proto (cutpoint_rdr); + struct casewriter *w = + autopaging_writer_create (proto); + struct casereader *r = casereader_clone (cutpoint_rdr); + struct ccase *cpc; + double prev_cp = SYSMIS; + + for ( ; (cpc = casereader_read (r) ); case_unref (cpc)) + { + struct ccase *new_case; + const double cp = case_data_idx (cpc, CUTPOINT)->f; + + assert (cp != SYSMIS); + + /* We don't want duplicates here */ + if ( cp == prev_cp ) + continue; + + new_case = case_clone (cpc); + + if ( pos_cond (result, cp)) + case_data_rw_idx (new_case, true_index)->f += weight; + else + case_data_rw_idx (new_case, false_index)->f += weight; + + prev_cp = cp; + + casewriter_write (w, new_case); + } + casereader_destroy (r); + + return casewriter_make_reader (w); +} + + + +static void output_roc (struct roc_state *rs, const struct cmd_roc *roc); + +/* + This function does 3 things: + + 1. Counts the number of cases which are equal to every other case in READER, + and those cases for which the relationship between it and every other case + satifies PRED (normally either > or <). VAR is variable defining a case's value + for this purpose. + + 2. Counts the number of true and false cases in reader, and populates + CUTPOINT_RDR accordingly. TRUE_INDEX and FALSE_INDEX are the indices + which receive these values. POS_COND is the condition defining true + and false. + + 3. CC is filled with the cumulative weight of all cases of READER. +*/ +static struct casereader * +process_group (const struct variable *var, struct casereader *reader, + bool (*pred) (double, double), + const struct dictionary *dict, + double *cc, + struct casereader **cutpoint_rdr, + bool (*pos_cond) (double, double), + int true_index, + int false_index) +{ + const struct variable *w = dict_get_weight (dict); + + struct casereader *r1 = + casereader_create_distinct (sort_execute_1var (reader, var), var, w); + + const int weight_idx = w ? var_get_case_index (w) : + caseproto_get_n_widths (casereader_get_proto (r1)) - 1; + + struct ccase *c1; + + struct casereader *rclone = casereader_clone (r1); + struct casewriter *wtr; + struct caseproto *proto = caseproto_create (); + + proto = caseproto_add_width (proto, 0); + proto = caseproto_add_width (proto, 0); + proto = caseproto_add_width (proto, 0); + + wtr = autopaging_writer_create (proto); + + *cc = 0; + + for ( ; (c1 = casereader_read (r1) ); case_unref (c1)) + { + struct ccase *new_case = case_create (proto); + struct ccase *c2; + struct casereader *r2 = casereader_clone (rclone); + + const double weight1 = case_data_idx (c1, weight_idx)->f; + const double d1 = case_data (c1, var)->f; + double n_eq = 0.0; + double n_pred = 0.0; + + *cutpoint_rdr = accumulate_counts (*cutpoint_rdr, d1, weight1, + pos_cond, + true_index, false_index); + + *cc += weight1; + + for ( ; (c2 = casereader_read (r2) ); case_unref (c2)) + { + const double d2 = case_data (c2, var)->f; + const double weight2 = case_data_idx (c2, weight_idx)->f; + + if ( d1 == d2 ) + { + n_eq += weight2; + continue; + } + else if ( pred (d2, d1)) + { + n_pred += weight2; + } + } + + case_data_rw_idx (new_case, VALUE)->f = d1; + case_data_rw_idx (new_case, N_EQ)->f = n_eq; + case_data_rw_idx (new_case, N_PRED)->f = n_pred; + + casewriter_write (wtr, new_case); + + casereader_destroy (r2); + } + + casereader_destroy (r1); + casereader_destroy (rclone); + + return casewriter_make_reader (wtr); +} + +/* Some more indeces into case data */ +#define N_POS_EQ 1 /* number of positive cases with values equal to n */ +#define N_POS_GT 2 /* number of postive cases with values greater than n */ +#define N_NEG_EQ 3 /* number of negative cases with values equal to n */ +#define N_NEG_LT 4 /* number of negative cases with values less than n */ + +static bool +gt (double d1, double d2) +{ + return d1 > d2; +} + + +static bool +ge (double d1, double d2) +{ + return d1 > d2; +} + +static bool +lt (double d1, double d2) +{ + return d1 < d2; +} + + +/* + Return a casereader with width 3, + populated with cases based upon READER. + The cases will have the values: + (N, number of cases equal to N, number of cases greater than N) + As a side effect, update RS->n1 with the number of positive cases. +*/ +static struct casereader * +process_positive_group (const struct variable *var, struct casereader *reader, + const struct dictionary *dict, + struct roc_state *rs) +{ + return process_group (var, reader, gt, dict, &rs->n1, + &rs->cutpoint_rdr, + ge, + TP, FN); +} + +/* + Return a casereader with width 3, + populated with cases based upon READER. + The cases will have the values: + (N, number of cases equal to N, number of cases less than N) + As a side effect, update RS->n2 with the number of negative cases. +*/ +static struct casereader * +process_negative_group (const struct variable *var, struct casereader *reader, + const struct dictionary *dict, + struct roc_state *rs) +{ + return process_group (var, reader, lt, dict, &rs->n2, + &rs->cutpoint_rdr, + lt, + TN, FP); +} + + + + +static void +append_cutpoint (struct casewriter *writer, double cutpoint) +{ + struct ccase *cc = case_create (casewriter_get_proto (writer)); + + case_data_rw_idx (cc, CUTPOINT)->f = cutpoint; + case_data_rw_idx (cc, TP)->f = 0; + case_data_rw_idx (cc, FN)->f = 0; + case_data_rw_idx (cc, TN)->f = 0; + case_data_rw_idx (cc, FP)->f = 0; + + casewriter_write (writer, cc); +} + + +/* + Create and initialise the rs[x].cutpoint_rdr casereaders. That is, the readers will + be created with width 5, ready to take the values (cutpoint, TP, FN, TN, FP), and the + reader will be populated with its final number of cases. + However on exit from this function, only CUTPOINT entries will be set to their final + value. The other entries will be initialised to zero. +*/ +static void +prepare_cutpoints (struct cmd_roc *roc, struct roc_state *rs, struct casereader *input) +{ + int i; + struct casereader *r = casereader_clone (input); + struct ccase *c; + struct caseproto *proto = caseproto_create (); + + struct subcase ordering; + subcase_init (&ordering, CUTPOINT, 0, SC_ASCEND); + + proto = caseproto_add_width (proto, 0); /* cutpoint */ + proto = caseproto_add_width (proto, 0); /* TP */ + proto = caseproto_add_width (proto, 0); /* FN */ + proto = caseproto_add_width (proto, 0); /* TN */ + proto = caseproto_add_width (proto, 0); /* FP */ + + for (i = 0 ; i < roc->n_vars; ++i) + { + rs[i].cutpoint_wtr = sort_create_writer (&ordering, proto); + rs[i].prev_result = SYSMIS; + rs[i].max = -DBL_MAX; + rs[i].min = DBL_MAX; + } + + for (; (c = casereader_read (r)) != NULL; case_unref (c)) + { + for (i = 0 ; i < roc->n_vars; ++i) + { + const union value *v = case_data (c, roc->vars[i]); + const double result = v->f; + + if ( mv_is_value_missing (var_get_missing_values (roc->vars[i]), v, roc->exclude)) + continue; + + minimize (&rs[i].min, result); + maximize (&rs[i].max, result); + + if ( rs[i].prev_result != SYSMIS && rs[i].prev_result != result ) + { + const double mean = (result + rs[i].prev_result ) / 2.0; + append_cutpoint (rs[i].cutpoint_wtr, mean); + } + + rs[i].prev_result = result; + } + } + casereader_destroy (r); + + + /* Append the min and max cutpoints */ + for (i = 0 ; i < roc->n_vars; ++i) + { + append_cutpoint (rs[i].cutpoint_wtr, rs[i].min - 1); + append_cutpoint (rs[i].cutpoint_wtr, rs[i].max + 1); + + rs[i].cutpoint_rdr = casewriter_make_reader (rs[i].cutpoint_wtr); + } +} + +static void +do_roc (struct cmd_roc *roc, struct casereader *reader, struct dictionary *dict) +{ + int i; + + struct roc_state *rs = xcalloc (roc->n_vars, sizeof *rs); + + struct casereader *negatives = NULL; + struct casereader *positives = NULL; + + struct caseproto *n_proto = caseproto_create (); + + struct subcase up_ordering; + struct subcase down_ordering; + + struct casewriter *neg_wtr = NULL; + + struct casereader *input = casereader_create_filter_missing (reader, + roc->vars, roc->n_vars, + roc->exclude, + NULL, + NULL); + + input = casereader_create_filter_missing (input, + &roc->state_var, 1, + roc->exclude, + NULL, + NULL); + + neg_wtr = autopaging_writer_create (casereader_get_proto (input)); + + prepare_cutpoints (roc, rs, input); + + + /* Separate the positive actual state cases from the negative ones */ + positives = + casereader_create_filter_func (input, + match_positives, + NULL, + roc, + neg_wtr); + + n_proto = caseproto_create (); + + n_proto = caseproto_add_width (n_proto, 0); + n_proto = caseproto_add_width (n_proto, 0); + n_proto = caseproto_add_width (n_proto, 0); + n_proto = caseproto_add_width (n_proto, 0); + n_proto = caseproto_add_width (n_proto, 0); + + subcase_init (&up_ordering, VALUE, 0, SC_ASCEND); + subcase_init (&down_ordering, VALUE, 0, SC_DESCEND); + + for (i = 0 ; i < roc->n_vars; ++i) + { + struct casewriter *w = NULL; + struct casereader *r = NULL; + + struct ccase *c; + + struct ccase *cpos; + struct casereader *n_neg ; + const struct variable *var = roc->vars[i]; + + struct casereader *neg ; + struct casereader *pos = casereader_clone (positives); + + + struct casereader *n_pos = + process_positive_group (var, pos, dict, &rs[i]); + + if ( negatives == NULL) + { + negatives = casewriter_make_reader (neg_wtr); + } + + neg = casereader_clone (negatives); + + n_neg = process_negative_group (var, neg, dict, &rs[i]); + + + /* Merge the n_pos and n_neg casereaders */ + w = sort_create_writer (&up_ordering, n_proto); + for ( ; (cpos = casereader_read (n_pos) ); case_unref (cpos)) + { + struct ccase *pos_case = case_create (n_proto); + struct ccase *cneg; + const double jpos = case_data_idx (cpos, VALUE)->f; + + while ((cneg = casereader_read (n_neg))) + { + struct ccase *nc = case_create (n_proto); + + const double jneg = case_data_idx (cneg, VALUE)->f; + + case_data_rw_idx (nc, VALUE)->f = jneg; + case_data_rw_idx (nc, N_POS_EQ)->f = 0; + + case_data_rw_idx (nc, N_POS_GT)->f = SYSMIS; + + *case_data_rw_idx (nc, N_NEG_EQ) = *case_data_idx (cneg, N_EQ); + *case_data_rw_idx (nc, N_NEG_LT) = *case_data_idx (cneg, N_PRED); + + casewriter_write (w, nc); + + case_unref (cneg); + if ( jneg > jpos) + break; + } + + case_data_rw_idx (pos_case, VALUE)->f = jpos; + *case_data_rw_idx (pos_case, N_POS_EQ) = *case_data_idx (cpos, N_EQ); + *case_data_rw_idx (pos_case, N_POS_GT) = *case_data_idx (cpos, N_PRED); + case_data_rw_idx (pos_case, N_NEG_EQ)->f = 0; + case_data_rw_idx (pos_case, N_NEG_LT)->f = SYSMIS; + + casewriter_write (w, pos_case); + } + +/* These aren't used anymore */ +#undef N_EQ +#undef N_PRED + + r = casewriter_make_reader (w); + + /* Propagate the N_POS_GT values from the positive cases + to the negative ones */ + { + double prev_pos_gt = rs[i].n1; + w = sort_create_writer (&down_ordering, n_proto); + + for ( ; (c = casereader_read (r) ); case_unref (c)) + { + double n_pos_gt = case_data_idx (c, N_POS_GT)->f; + struct ccase *nc = case_clone (c); + + if ( n_pos_gt == SYSMIS) + { + n_pos_gt = prev_pos_gt; + case_data_rw_idx (nc, N_POS_GT)->f = n_pos_gt; + } + + casewriter_write (w, nc); + prev_pos_gt = n_pos_gt; + } + + r = casewriter_make_reader (w); + } + + /* Propagate the N_NEG_LT values from the negative cases + to the positive ones */ + { + double prev_neg_lt = rs[i].n2; + w = sort_create_writer (&up_ordering, n_proto); + + for ( ; (c = casereader_read (r) ); case_unref (c)) + { + double n_neg_lt = case_data_idx (c, N_NEG_LT)->f; + struct ccase *nc = case_clone (c); + + if ( n_neg_lt == SYSMIS) + { + n_neg_lt = prev_neg_lt; + case_data_rw_idx (nc, N_NEG_LT)->f = n_neg_lt; + } + + casewriter_write (w, nc); + prev_neg_lt = n_neg_lt; + } + + r = casewriter_make_reader (w); + } + + { + struct ccase *prev_case = NULL; + for ( ; (c = casereader_read (r) ); case_unref (c)) + { + const struct ccase *next_case = casereader_peek (r, 0); + + const double j = case_data_idx (c, VALUE)->f; + double n_pos_eq = case_data_idx (c, N_POS_EQ)->f; + double n_pos_gt = case_data_idx (c, N_POS_GT)->f; + double n_neg_eq = case_data_idx (c, N_NEG_EQ)->f; + double n_neg_lt = case_data_idx (c, N_NEG_LT)->f; + + if ( prev_case && j == case_data_idx (prev_case, VALUE)->f) + { + if ( 0 == case_data_idx (c, N_POS_EQ)->f) + { + n_pos_eq = case_data_idx (prev_case, N_POS_EQ)->f; + n_pos_gt = case_data_idx (prev_case, N_POS_GT)->f; + } + + if ( 0 == case_data_idx (c, N_NEG_EQ)->f) + { + n_neg_eq = case_data_idx (prev_case, N_NEG_EQ)->f; + n_neg_lt = case_data_idx (prev_case, N_NEG_LT)->f; + } + } + + if ( NULL == next_case || j != case_data_idx (next_case, VALUE)->f) + { + rs[i].auc += n_pos_gt * n_neg_eq + (n_pos_eq * n_neg_eq) / 2.0; + + rs[i].q1hat += + n_neg_eq * ( pow2 (n_pos_gt) + n_pos_gt * n_pos_eq + pow2 (n_pos_eq) / 3.0); + rs[i].q2hat += + n_pos_eq * ( pow2 (n_neg_lt) + n_neg_lt * n_neg_eq + pow2 (n_neg_eq) / 3.0); + + } + + case_unref (prev_case); + prev_case = case_clone (c); + } + + rs[i].auc /= rs[i].n1 * rs[i].n2; + if ( roc->invert ) + rs[i].auc = 1 - rs[i].auc; + + if ( roc->bi_neg_exp ) + { + rs[i].q1hat = rs[i].auc / ( 2 - rs[i].auc); + rs[i].q2hat = 2 * pow2 (rs[i].auc) / ( 1 + rs[i].auc); + } + else + { + rs[i].q1hat /= rs[i].n2 * pow2 (rs[i].n1); + rs[i].q2hat /= rs[i].n1 * pow2 (rs[i].n2); + } + } + } + + casereader_destroy (positives); + casereader_destroy (negatives); + + output_roc (rs, roc); + + free (rs); +} + +static void +show_auc (struct roc_state *rs, const struct cmd_roc *roc) +{ + int i; + const int n_fields = roc->print_se ? 5 : 1; + const int n_cols = roc->n_vars > 1 ? n_fields + 1: n_fields; + const int n_rows = 2 + roc->n_vars; + struct tab_table *tbl = tab_create (n_cols, n_rows, 0); + + if ( roc->n_vars > 1) + tab_title (tbl, _("Area Under the Curve")); + else + tab_title (tbl, _("Area Under the Curve (%s)"), var_to_string (roc->vars[0])); + + tab_headers (tbl, n_cols - n_fields, 0, 1, 0); + + tab_dim (tbl, tab_natural_dimensions, NULL); + + tab_text (tbl, n_cols - n_fields, 1, TAT_TITLE, _("Area")); + + tab_hline (tbl, TAL_2, 0, n_cols - 1, 2); + + tab_box (tbl, + TAL_2, TAL_2, + -1, TAL_1, + 0, 0, + n_cols - 1, + n_rows - 1); + + if ( roc->print_se ) + { + tab_text (tbl, n_cols - 4, 1, TAT_TITLE, _("Std. Error")); + tab_text (tbl, n_cols - 3, 1, TAT_TITLE, _("Asymptotic Sig.")); + + tab_text (tbl, n_cols - 2, 1, TAT_TITLE, _("Lower Bound")); + tab_text (tbl, n_cols - 1, 1, TAT_TITLE, _("Upper Bound")); + + tab_joint_text (tbl, n_cols - 2, 0, 4, 0, + TAT_TITLE | TAB_CENTER | TAT_PRINTF, + _("Asymp. %g%% Confidence Interval"), roc->ci); + tab_vline (tbl, 0, n_cols - 1, 0, 0); + tab_hline (tbl, TAL_1, n_cols - 2, n_cols - 1, 1); + } + + if ( roc->n_vars > 1) + tab_text (tbl, 0, 1, TAT_TITLE, _("Variable under test")); + + if ( roc->n_vars > 1) + tab_vline (tbl, TAL_2, 1, 0, n_rows - 1); + + + for ( i = 0 ; i < roc->n_vars ; ++i ) + { + tab_text (tbl, 0, 2 + i, TAT_TITLE, var_to_string (roc->vars[i])); + + tab_double (tbl, n_cols - n_fields, 2 + i, 0, rs[i].auc, NULL); + + if ( roc->print_se ) + { + double se ; + const double sd_0_5 = sqrt ((rs[i].n1 + rs[i].n2 + 1) / + (12 * rs[i].n1 * rs[i].n2)); + double ci ; + double yy ; + + se = rs[i].auc * (1 - rs[i].auc) + (rs[i].n1 - 1) * (rs[i].q1hat - pow2 (rs[i].auc)) + + (rs[i].n2 - 1) * (rs[i].q2hat - pow2 (rs[i].auc)); + + se /= rs[i].n1 * rs[i].n2; + + se = sqrt (se); + + tab_double (tbl, n_cols - 4, 2 + i, 0, + se, + NULL); + + ci = 1 - roc->ci / 100.0; + yy = gsl_cdf_gaussian_Qinv (ci, se) ; + + tab_double (tbl, n_cols - 2, 2 + i, 0, + rs[i].auc - yy, + NULL); + + tab_double (tbl, n_cols - 1, 2 + i, 0, + rs[i].auc + yy, + NULL); + + tab_double (tbl, n_cols - 3, 2 + i, 0, + 2.0 * gsl_cdf_ugaussian_Q (fabs ((rs[i].auc - 0.5 ) / sd_0_5)), + NULL); + } + } + + tab_submit (tbl); +} + + +static void +show_summary (const struct cmd_roc *roc) +{ + const int n_cols = 3; + const int n_rows = 4; + struct tab_table *tbl = tab_create (n_cols, n_rows, 0); + + tab_title (tbl, _("Case Summary")); + + tab_headers (tbl, 1, 0, 2, 0); + + tab_dim (tbl, tab_natural_dimensions, NULL); + + tab_box (tbl, + TAL_2, TAL_2, + -1, -1, + 0, 0, + n_cols - 1, + n_rows - 1); + + tab_hline (tbl, TAL_2, 0, n_cols - 1, 2); + tab_vline (tbl, TAL_2, 1, 0, n_rows - 1); + + + tab_hline (tbl, TAL_2, 1, n_cols - 1, 1); + tab_vline (tbl, TAL_1, 2, 1, n_rows - 1); + + + tab_text (tbl, 0, 1, TAT_TITLE | TAB_LEFT, var_to_string (roc->state_var)); + tab_text (tbl, 1, 1, TAT_TITLE, _("Unweighted")); + tab_text (tbl, 2, 1, TAT_TITLE, _("Weighted")); + + tab_joint_text (tbl, 1, 0, 2, 0, + TAT_TITLE | TAB_CENTER, + _("Valid N (listwise)")); + + + tab_text (tbl, 0, 2, TAB_LEFT, _("Positive")); + tab_text (tbl, 0, 3, TAB_LEFT, _("Negative")); + + + tab_double (tbl, 1, 2, 0, roc->pos, &F_8_0); + tab_double (tbl, 1, 3, 0, roc->neg, &F_8_0); + + tab_double (tbl, 2, 2, 0, roc->pos_weighted, 0); + tab_double (tbl, 2, 3, 0, roc->neg_weighted, 0); + + tab_submit (tbl); +} + + +static void +show_coords (struct roc_state *rs, const struct cmd_roc *roc) +{ + int x = 1; + int i; + const int n_cols = roc->n_vars > 1 ? 4 : 3; + int n_rows = 1; + struct tab_table *tbl ; + + for (i = 0; i < roc->n_vars; ++i) + n_rows += casereader_count_cases (rs[i].cutpoint_rdr); + + tbl = tab_create (n_cols, n_rows, 0); + + if ( roc->n_vars > 1) + tab_title (tbl, _("Coordinates of the Curve")); + else + tab_title (tbl, _("Coordinates of the Curve (%s)"), var_to_string (roc->vars[0])); + + + tab_headers (tbl, 1, 0, 1, 0); + + tab_dim (tbl, tab_natural_dimensions, NULL); + + tab_hline (tbl, TAL_2, 0, n_cols - 1, 1); + + if ( roc->n_vars > 1) + tab_text (tbl, 0, 0, TAT_TITLE, _("Test variable")); + + tab_text (tbl, n_cols - 3, 0, TAT_TITLE, _("Positive if greater than or equal to")); + tab_text (tbl, n_cols - 2, 0, TAT_TITLE, _("Sensitivity")); + tab_text (tbl, n_cols - 1, 0, TAT_TITLE, _("1 - Specificity")); + + tab_box (tbl, + TAL_2, TAL_2, + -1, TAL_1, + 0, 0, + n_cols - 1, + n_rows - 1); + + if ( roc->n_vars > 1) + tab_vline (tbl, TAL_2, 1, 0, n_rows - 1); + + for (i = 0; i < roc->n_vars; ++i) + { + struct ccase *cc; + struct casereader *r = casereader_clone (rs[i].cutpoint_rdr); + + if ( roc->n_vars > 1) + tab_text (tbl, 0, x, TAT_TITLE, var_to_string (roc->vars[i])); + + if ( i > 0) + tab_hline (tbl, TAL_1, 0, n_cols - 1, x); + + + for (; (cc = casereader_read (r)) != NULL; + case_unref (cc), x++) + { + const double se = case_data_idx (cc, TP)->f / + ( + case_data_idx (cc, TP)->f + + + case_data_idx (cc, FN)->f + ); + + const double sp = case_data_idx (cc, TN)->f / + ( + case_data_idx (cc, TN)->f + + + case_data_idx (cc, FP)->f + ); + + tab_double (tbl, n_cols - 3, x, 0, case_data_idx (cc, CUTPOINT)->f, + var_get_print_format (roc->vars[i])); + + tab_double (tbl, n_cols - 2, x, 0, se, NULL); + tab_double (tbl, n_cols - 1, x, 0, 1 - sp, NULL); + } + + casereader_destroy (r); + } + + tab_submit (tbl); +} + + +static void +draw_roc (struct roc_state *rs, const struct cmd_roc *roc) +{ + int i; + + struct chart *roc_chart = chart_create (); + + chart_write_title (roc_chart, _("ROC Curve")); + chart_write_xlabel (roc_chart, _("1 - Specificity")); + chart_write_ylabel (roc_chart, _("Sensitivity")); + + chart_write_xscale (roc_chart, 0, 1, 5); + chart_write_yscale (roc_chart, 0, 1, 5); + + if ( roc->reference ) + { + chart_line (roc_chart, 1.0, 0, + 0.0, 1.0, + CHART_DIM_X); + } + + for (i = 0; i < roc->n_vars; ++i) + { + struct ccase *cc; + struct casereader *r = casereader_clone (rs[i].cutpoint_rdr); + + chart_vector_start (roc_chart, var_get_name (roc->vars[i])); + for (; (cc = casereader_read (r)) != NULL; + case_unref (cc)) + { + double se = case_data_idx (cc, TP)->f; + double sp = case_data_idx (cc, TN)->f; + + se /= case_data_idx (cc, FN)->f + + case_data_idx (cc, TP)->f ; + + sp /= case_data_idx (cc, TN)->f + + case_data_idx (cc, FP)->f ; + + chart_vector (roc_chart, 1 - sp, se); + } + chart_vector_end (roc_chart); + casereader_destroy (r); + } + + chart_write_legend (roc_chart); + + chart_submit (roc_chart); +} + + +static void +output_roc (struct roc_state *rs, const struct cmd_roc *roc) +{ + show_summary (roc); + + if ( roc->curve ) + draw_roc (rs, roc); + + show_auc (rs, roc); + + + if ( roc->print_coords ) + show_coords (rs, roc); +} + diff --git a/src/output/chart.c b/src/output/chart.c index 49819e9f..e324901c 100644 --- a/src/output/chart.c +++ b/src/output/chart.c @@ -83,6 +83,9 @@ chart_create(void) chart->legend_left = 810; chart->legend_right = 1000; chart->font_size = 0; + chart->in_path = false; + chart->dataset = NULL; + chart->n_datasets = 0; strcpy(chart->fill_colour,"red"); /* Get default font size */ @@ -100,6 +103,7 @@ chart_create(void) void chart_submit(struct chart *chart) { + int i; struct som_entity s; struct outp_driver *d; @@ -124,6 +128,11 @@ chart_submit(struct chart *chart) d = outp_drivers (NULL); d->class->finalise_chart(d, chart); + + for (i = 0 ; i < chart->n_datasets; ++i) + free (chart->dataset[i]); + free (chart->dataset); + free(chart); } diff --git a/src/output/chart.h b/src/output/chart.h index 9586eb1e..4da12ecf 100644 --- a/src/output/chart.h +++ b/src/output/chart.h @@ -64,6 +64,8 @@ struct chart { int legend_left ; int legend_right ; + const char **dataset; + int n_datasets; /* Default font size for the plot (if zero, then use plotter default) */ @@ -78,6 +80,7 @@ struct chart { double x_max; double y_min; double y_max; + bool in_path; }; diff --git a/src/output/charts/cartesian.c b/src/output/charts/cartesian.c index 75c49aa8..cb2f346b 100644 --- a/src/output/charts/cartesian.c +++ b/src/output/charts/cartesian.c @@ -27,26 +27,26 @@ #include - -struct dataset +/* Start a new vector called NAME */ +void +chart_vector_start (struct chart *ch, const char *name) { - int n_data; - const char *label; -}; + if ( ! ch ) + return ; + pl_savestate_r (ch->lp); -#define DATASETS 2 + pl_colorname_r (ch->lp, data_colour [ch->n_datasets % N_CHART_COLOURS]); -static const struct dataset dataset[DATASETS] = - { - { 13, "male"}, - { 11, "female"}, - }; + ch->n_datasets++; + ch->dataset = xrealloc (ch->dataset, ch->n_datasets * sizeof (*ch->dataset)); + ch->dataset[ch->n_datasets - 1] = strdup (name); +} /* Plot a data point */ void -chart_datum(struct chart *ch, int dataset UNUSED, double x, double y) +chart_datum (struct chart *ch, int dataset UNUSED, double x, double y) { if ( ! ch ) return ; @@ -66,13 +66,48 @@ chart_datum(struct chart *ch, int dataset UNUSED, double x, double y) } } +void +chart_vector_end (struct chart *ch) +{ + pl_endpath_r (ch->lp); + pl_colorname_r (ch->lp, "black"); + ch->in_path = false; + pl_restorestate_r (ch->lp); +} + +/* Plot a data point */ +void +chart_vector (struct chart *ch, double x, double y) +{ + if ( ! ch ) + return ; + + { + const double x_pos = + (x - ch->x_min) * ch->abscissa_scale + ch->data_left ; + + const double y_pos = + (y - ch->y_min) * ch->ordinate_scale + ch->data_bottom ; + + if ( ch->in_path) + pl_fcont_r (ch->lp, x_pos, y_pos); + else + { + pl_fmove_r (ch->lp, x_pos, y_pos); + ch->in_path = true; + } + } +} + + + /* Draw a line with slope SLOPE and intercept INTERCEPT. between the points limit1 and limit2. If lim_dim is CHART_DIM_Y then the limit{1,2} are on the y axis otherwise the x axis */ void -chart_line(struct chart *ch, double slope, double intercept, +chart_line (struct chart *ch, double slope, double intercept, double limit1, double limit2, enum CHART_DIM lim_dim) { double x1, y1; @@ -107,5 +142,4 @@ chart_line(struct chart *ch, double slope, double intercept, pl_fline_r(ch->lp, x1, y1, x2, y2); pl_restorestate_r(ch->lp); - } diff --git a/src/output/charts/cartesian.h b/src/output/charts/cartesian.h index 40e3f714..0874b9cc 100644 --- a/src/output/charts/cartesian.h +++ b/src/output/charts/cartesian.h @@ -27,16 +27,19 @@ enum CHART_DIM }; +void chart_vector_start (struct chart *ch, const char *name); +void chart_vector (struct chart *ch, double x, double y); +void chart_vector_end (struct chart *ch); /* Plot a data point */ -void chart_datum(struct chart *ch, int dataset UNUSED, double x, double y); +void chart_datum (struct chart *ch, int dataset UNUSED, double x, double y); /* Draw a line with slope SLOPE and intercept INTERCEPT. between the points limit1 and limit2. If lim_dim is CHART_DIM_Y then the limit{1,2} are on the y axis otherwise the x axis */ -void chart_line(struct chart *ch, double slope, double intercept, +void chart_line (struct chart *ch, double slope, double intercept, double limit1, double limit2, enum CHART_DIM lim_dim); diff --git a/src/output/charts/plot-chart.c b/src/output/charts/plot-chart.c index 3b4f1b37..5641db12 100644 --- a/src/output/charts/plot-chart.c +++ b/src/output/charts/plot-chart.c @@ -1,5 +1,5 @@ /* PSPP - a program for statistical analysis. - Copyright (C) 2004 Free Software Foundation, Inc. + Copyright (C) 2004, 2009 Free Software Foundation, Inc. This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -210,3 +210,47 @@ chart_write_ylabel(struct chart *ch, const char *label) pl_restorestate_r(ch->lp); } + + +void +chart_write_legend (struct chart *ch) +{ + int i; + const int vstep = ch->font_size * 2; + const int xpad = 10; + const int ypad = 10; + const int swatch = 20; + const int legend_top = ch->data_top; + const int legend_bottom = legend_top - + (vstep * ch->n_datasets + 2 * ypad ); + + if ( ! ch ) + return ; + + pl_savestate_r (ch->lp); + + pl_box_r (ch->lp, ch->legend_left, legend_top, + ch->legend_right - xpad, legend_bottom); + + for (i = 0 ; i < ch->n_datasets ; ++i ) + { + const int ypos = vstep * (i + 1); + const int xpos = ch->legend_left + xpad; + pl_move_r (ch->lp, xpos, legend_top - ypos); + + pl_savestate_r(ch->lp); + pl_fillcolorname_r (ch->lp, data_colour [ i % N_CHART_COLOURS]); + + pl_pentype_r (ch->lp, 1); + pl_filltype_r (ch->lp, 1); + pl_boxrel_r (ch->lp, 0, 0, swatch, swatch); + + + pl_moverel_r (ch->lp, swatch, 0); + pl_alabel_r (ch->lp, 0, 0, ch->dataset[i]); + + pl_restorestate_r (ch->lp); + } + + pl_restorestate_r (ch->lp); +} diff --git a/src/output/charts/plot-chart.h b/src/output/charts/plot-chart.h index 4a1dc105..f4cc5bbf 100644 --- a/src/output/charts/plot-chart.h +++ b/src/output/charts/plot-chart.h @@ -70,4 +70,6 @@ void chart_write_xlabel(struct chart *ch, const char *label) ; /* Write the ordinate label */ void chart_write_ylabel(struct chart *ch, const char *label); +void chart_write_legend (struct chart *ch); + #endif diff --git a/tests/automake.mk b/tests/automake.mk index 7bfd77bf..aba1a62c 100644 --- a/tests/automake.mk +++ b/tests/automake.mk @@ -56,6 +56,8 @@ dist_TESTS = \ tests/command/regression.sh \ tests/command/regression-qr.sh \ tests/command/reliability.sh \ + tests/command/roc.sh \ + tests/command/roc2.sh \ tests/command/sample.sh \ tests/command/sort.sh \ tests/command/sysfiles.sh \ diff --git a/tests/command/roc.sh b/tests/command/roc.sh new file mode 100755 index 00000000..0c24f0c5 --- /dev/null +++ b/tests/command/roc.sh @@ -0,0 +1,173 @@ +#!/bin/sh + +# This program tests the ROC command. + +TEMPDIR=/tmp/pspp-tst-$$ +TESTFILE=$TEMPDIR/`basename $0`.sps + +# ensure that top_builddir are absolute +if [ -z "$top_builddir" ] ; then top_builddir=. ; fi +if [ -z "$top_srcdir" ] ; then top_srcdir=. ; fi +top_builddir=`cd $top_builddir; pwd` +PSPP=$top_builddir/src/ui/terminal/pspp + +# ensure that top_srcdir is absolute +top_srcdir=`cd $top_srcdir; pwd` + +STAT_CONFIG_PATH=$top_srcdir/config +export STAT_CONFIG_PATH + +LANG=C +export LANG + +cleanup() +{ + if [ x"$PSPP_TEST_NO_CLEANUP" != x ] ; then + echo "NOT cleaning $TEMPDIR" + return ; + fi + cd / + rm -rf $TEMPDIR +} + + +fail() +{ + echo $activity + echo FAILED + cleanup; + exit 1; +} + + +no_result() +{ + echo $activity + echo NO RESULT; + cleanup; + exit 2; +} + +pass() +{ + cleanup; + exit 0; +} + +mkdir -p $TEMPDIR + +cd $TEMPDIR + +activity="create program" +cat > $TESTFILE < $TESTFILE <