realset_smob.cpp

非常著名的曲线拟合程序

// This is -*- C++ -*-
// $Id: RealSet_smob.cpp,v 1.13 1999/07/15 19:01:54 trow Exp $

/* 
 * RealSet_smob.cpp
 *
 * Copyright (C) 1999 EMC Capital Management
 *
 * Developed by Jon Trowbridge & Havoc Pennington
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library 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
 * Library General Public License for more details.
 *
 * You should have received a copy of the GNU Library General Public
 * License along with this library; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
 * 02111-1307, USA.
 */

#include <RealSet.h>
#include <descriptive.h>
#include "guilegoose.h"
#include "snarf.h"

// This is used to decide whether to delete the realset
//  when we free the smob
template <class T>
class OwnerPtr 
{
public:
  OwnerPtr(T* t) : ptr_(t), owned_(true) {}
  ~OwnerPtr() {}

  T* ptr() { return ptr_; }
  
  void disown() { owned_ = false; }

  bool owned() { return owned_; }

private:
  T* const ptr_;
  bool owned_;
};


static long realset_type_tag;

#define SCM_TO_OWNERREALSET(obj) (reinterpret_cast<OwnerPtr<RealSet>*>(SCM_CDR(obj)))
#define SCM_TO_REALSET(obj) (SCM_TO_OWNERREALSET(obj)->ptr())
#define REALSET_P(value)    (SCM_NIMP (value) && SCM_CAR(value) == realset_type_tag)

bool     
scm_realsetp (SCM obj)
{
  return REALSET_P(obj);
}

RealSet* 
scm2realset (SCM obj)
{
  if (!REALSET_P(obj)) return 0;
  return SCM_TO_REALSET(obj);
}

SCM      
realset2scm (RealSet* rs)
{
  SCM_DEFER_INTS;
  OwnerPtr<RealSet>* op = new OwnerPtr<RealSet>(rs);
  scm_done_malloc(sizeof(OwnerPtr<RealSet>));
  SCM smob;
  SCM_NEWCELL (smob);
  SCM_SETCDR (smob, op);
  SCM_SETCAR (smob, realset_type_tag);
  SCM_ALLOW_INTS;
  return smob;
}

SCM 
realset2scm_handle(RealSet* rs)
{
  SCM obj = realset2scm(rs);
  SCM_TO_OWNERREALSET(obj)->disown(); // so we won't destroy it.
  return obj;
}

// This function should mark any SCM objects we reference,
//  at this time none. The return value is automatically marked 
//  by Guile, this is some kind of optimization; so we can
//  return one of the objects.
static SCM
mark_realset (SCM obj)
{
  //  scm_gc_mark (image->name);
  //  return image->update_func;
  return SCM_BOOL_F; // not important
}

static scm_sizet
free_realset (SCM obj)
{
  OwnerPtr<RealSet>* op = SCM_TO_OWNERREALSET(obj);

  // we return the number of bytes freed, for statistics;
  //  hopefully it doesn't have to be exact...
  scm_sizet size = sizeof(OwnerPtr<RealSet>); 

  SCM_DEFER_INTS;

  if (op->owned())
    {
      RealSet* rs = op->ptr();
      delete op;
      delete rs;
      size += sizeof (RealSet);
    }

  SCM_ALLOW_INTS;

  return size;
}

static int
print_realset (SCM obj, SCM port, scm_print_state *pstate)
{
  RealSet* rs = SCM_TO_REALSET(obj);

  string rep = "<";

  if (rs->label().size()) {
    rep += "RealSet `";
    rep += rs->label();
    rep += "'";
  } else {
    rep += "Unnamed RealSet";
  }
  rep += " (";
  if (rs->size()) {
    char buf[128];
    sprintf(buf, "size=%u", rs->size());
    rep += buf;
  } else {
    rep += "empty";
  }
  rep += ")>";

  scm_puts (const_cast<char*>(rep.c_str()), port);

  /* non-zero means success */
  return 1;
}

static scm_smobfuns realset_funcs = {
  mark_realset, 
  free_realset, 
  print_realset, 
  0 // means we can never be equal? 
};

//////////////////////////////////////////////////////////////////////////////

GOOSE_PROC(realsetp, "realset?",
	   1,0,0, (SCM x))
{
  return gh_bool2scm(REALSET_P(x));
}

GOOSE_PROC(make_realset, "make-realset",
	   0,0,1, (SCM scm_l))
{   
  SCM_DEFER_INTS;

  RealSet* rs = new RealSet;

  // Note that we would NOT do this if merely creating a handle.
  scm_done_malloc(sizeof(RealSet));

  SCM_ALLOW_INTS;

  if (! SCM_UNBNDP(scm_l)) {
    while(!gh_null_p(scm_l)) {
      SCM scm_x = gh_car(scm_l);
      SCM_ASSERT(gh_number_p(scm_x), scm_x,
		 "Non-numeric values can't be added to a realset",
		 str_make_realset);
      double x = gh_scm2double(scm_x);
      rs->add(x);
      scm_l = gh_cdr(scm_l);
    }
  }

  return realset2scm(rs);
}

GOOSE_PROC(rs_label, "realset-label",
	   1,0,0, (SCM scm_rs))
{
  SCM_ASSERT(REALSET_P(scm_rs), scm_rs, SCM_ARG1, str_rs_label);
  RealSet* rs = scm2realset(scm_rs);
  return gh_str02scm((char*)rs->label().c_str());
}

GOOSE_PROC(rs_set_label, "realset-set-label!",
	   2,0,0, (SCM scm_rs, SCM scm_lab))
{
  SCM_ASSERT(REALSET_P(scm_rs), scm_rs, SCM_ARG1, str_rs_set_label);
  SCM_ASSERT(gh_string_p(scm_lab), scm_lab, SCM_ARG2, str_rs_set_label);
  int len;
  char* foo = gh_scm2newstr(scm_lab, &len);
  string lab(foo);
  free(foo);
  scm2realset(scm_rs)->set_label(lab);
  return scm_rs;
}

GOOSE_PROC(realset_clone, "realset-clone",
	   1,0,0, (SCM scm_rs))
{
  SCM_ASSERT(REALSET_P(scm_rs), scm_rs, SCM_ARG1, str_realset_clone);
  RealSet* rs = new RealSet(*scm2realset(scm_rs));
  return realset2scm(rs);
}

GOOSE_PROC(realset_add, "realset-add!",
	   1,0,1, (SCM scm_rs, SCM scm_l))
{
  SCM_ASSERT(REALSET_P(scm_rs), scm_rs, SCM_ARG1, str_realset_add);
  
  RealSet* rs = scm2realset(scm_rs);

  while(!gh_null_p(scm_l)) {
    SCM scm_x = gh_car(scm_l);
    SCM_ASSERT(gh_number_p(scm_x), scm_x,
	       "Non-numeric arguments are not allowed in realset-add!",
	       str_realset_add);
    double x = gh_scm2double(scm_x);
    rs->add(x);
    scm_l = gh_cdr(scm_l);
  }
  
  return scm_rs;
}

GOOSE_PROC(realset_append, "realset-append!",
	   2,0,0, (SCM scm_rs, SCM scm_app))
{
  SCM_ASSERT(REALSET_P(scm_rs), scm_rs, SCM_ARG1, str_realset_append);
  SCM_ASSERT(REALSET_P(scm_app), scm_app, SCM_ARG2, str_realset_append);

  RealSet* rs = scm2realset(scm_rs);
  RealSet* app = scm2realset(scm_app);
  rs->append(*app);

  return scm_rs;
}

GOOSE_PROC(realset_size, "realset-size",
	   1,0,0, (SCM scm_rs))
{
  SCM_ASSERT(REALSET_P(scm_rs), scm_rs, SCM_ARG1, str_realset_size);
  RealSet* rs = scm2realset(scm_rs); 
  int n = rs->size();
  return gh_int2scm(n);
}

GOOSE_PROC(realset_ref, "realset-ref",
	   2,0,0, (SCM scm_rs, SCM scm_i))
{
  SCM_ASSERT(REALSET_P(scm_rs), scm_rs, SCM_ARG1, str_realset_ref);
  SCM_ASSERT(gh_exact_p(scm_i), scm_i, SCM_ARG2, str_realset_ref);
  RealSet* rs = scm2realset(scm_rs); 
  int i = gh_scm2int(scm_i);
  SCM_ASSERT(i>=0 && i < (int)(rs->size()), scm_i, SCM_OUTOFRANGE,
	     str_realset_ref);
  double x = rs->data(i);
  return gh_double2scm(x);
}

GOOSE_PROC(realset_sorted_ref, "realset-sorted-ref",
	   2,0,0, (SCM scm_rs, SCM scm_i))
{
  SCM_ASSERT(REALSET_P(scm_rs), scm_rs, SCM_ARG1, str_realset_sorted_ref);
  SCM_ASSERT(gh_exact_p(scm_i), scm_i, SCM_ARG2, str_realset_sorted_ref);
  RealSet* rs = scm2realset(scm_rs); 
  int i = gh_scm2int(scm_i);
  SCM_ASSERT(i>=0 || i < (int)(rs->size()), scm_i, SCM_OUTOFRANGE,
	     str_realset_sorted_ref);
  double x = rs->sorted_data(i);
  return gh_double2scm(x);
}

// Behold the perversity of macros within macros!
// Normally I'm opposed to this sort of thing, but in this case...
#define RSFUNC(fn, str) \
GOOSE_PROC(rsfunc_##fn, str, 1,0,0, (SCM scm_rs)) \
{ \
  SCM_ASSERT(REALSET_P(scm_rs), scm_rs, SCM_ARG1, str_rsfunc_##fn); \
  RealSet* rs = scm2realset(scm_rs); \
  double x=0; \
  try { x = goose_##fn(*rs); } GOOSE_CATCH(str_rsfunc_##fn);\
  return gh_double2scm(x); \
}

RSFUNC(min, "realset-min")
RSFUNC(max, "realset-max")
RSFUNC(range, "realset-range")
RSFUNC(sum, "realset-sum")
RSFUNC(mean, "realset-mean")
RSFUNC(var, "realset-var")
RSFUNC(sdev, "realset-sdev")
RSFUNC(sdevs, "realset-sdevs")

GOOSE_PROC(realset_percentile, "realset-percentile",
	   2,0,0, (SCM scm_rs, SCM scm_p))
{
  SCM_ASSERT(REALSET_P(scm_rs), scm_rs, SCM_ARG1, str_realset_percentile);
  SCM_ASSERT(gh_number_p(scm_p), scm_p, SCM_ARG2, str_realset_percentile);

  double x=0;
  try {
    RealSet* rs = scm2realset(scm_rs);
    double p = gh_scm2double(scm_p);
    x = goose_percentile(*rs, p);
  } GOOSE_CATCH(str_realset_percentile);

  return gh_double2scm(x);
}

RSFUNC(med, "realset-med")
RSFUNC(q1, "realset-q1")
RSFUNC(q3, "realset-q3")
RSFUNC(iqr, "realset-iqr")

GOOSE_PROC(realset_decile, "realset-decile",
	   2,0,0, (SCM scm_rs, SCM scm_i))
{
  SCM_ASSERT(REALSET_P(scm_rs), scm_rs, SCM_ARG1, str_realset_decile);
  SCM_ASSERT(gh_exact_p(scm_i), scm_i, SCM_ARG2, str_realset_decile);

  double x = 0;
  try {
    RealSet* rs = scm2realset(scm_rs);
    int i = gh_scm2int(scm_i);
    x = goose_decile(*rs, i);
  } GOOSE_CATCH(str_realset_decile);

  return gh_double2scm(x);
}

GOOSE_PROC(mean_trimmed, "realset-mean-trimmed",
	   2,1,0, (SCM scm_rs, SCM scm_l, SCM scm_r))
{
  SCM_ASSERT(REALSET_P(scm_rs), scm_rs, SCM_ARG1, str_mean_trimmed);
  SCM_ASSERT(gh_exact_p(scm_l), scm_l, SCM_ARG2, str_mean_trimmed);
  SCM_ASSERT(SCM_UNBNDP(scm_r) || gh_exact_p(scm_r),
	     scm_r, SCM_ARG2, str_mean_trimmed);
  double x = 0;
  try {
    RealSet* rs = scm2realset(scm_rs);
    int l = gh_scm2int(scm_l);
    int r = SCM_UNBNDP(scm_r) ? l : gh_scm2int(scm_r);
    x = goose_mean_trimmed(*rs,l,r);
  } GOOSE_CATCH(str_mean_trimmed);

  return gh_double2scm(x); 
}

GOOSE_PROC(mean_winsorized, "realset-mean-winsorized",
	   2,1,0, (SCM scm_rs, SCM scm_l, SCM scm_r))
{
  SCM_ASSERT(REALSET_P(scm_rs), scm_rs, SCM_ARG1, str_mean_winsorized);
  SCM_ASSERT(gh_exact_p(scm_l), scm_l, SCM_ARG2, str_mean_winsorized);
  SCM_ASSERT(SCM_UNBNDP(scm_r) || gh_exact_p(scm_r),
	     scm_r, SCM_ARG2, str_mean_winsorized);
  double x = 0;
  try {
    RealSet* rs = scm2realset(scm_rs);
    int l = gh_scm2int(scm_l);
    int r = SCM_UNBNDP(scm_r) ? l : gh_scm2int(scm_r);
    x = goose_mean_winsorized(*rs,l,r);
  } GOOSE_CATCH(str_mean_winsorized);

  return gh_double2scm(x); 
}

GOOSE_PROC(moment, "realset-moment",
	   2,1,0, (SCM scm_rs, SCM scm_k, SCM scm_x))
{
  SCM_ASSERT(REALSET_P(scm_rs), scm_rs, SCM_ARG1, str_moment);
  SCM_ASSERT(gh_exact_p(scm_k), scm_k, SCM_ARG2, str_moment);
  SCM_ASSERT(SCM_UNBNDP(scm_x) || gh_number_p(scm_x), scm_x, SCM_ARG3,
	     str_moment);
  double m = 0;
  try {
    RealSet* rs = scm2realset(scm_rs);
    int k = gh_scm2int(scm_k);
    double x = SCM_UNBNDP(scm_x) ? rs->mean() : gh_scm2double(scm_x);
    m = goose_moment(*rs,k,x);
  } GOOSE_CATCH(str_moment);

  return gh_double2scm(m);
}

RSFUNC(gmean, "realset-gmean")
RSFUNC(hmean, "realset-hmean")
RSFUNC(rms, "realset-rms")
RSFUNC(meandev, "realset-meandev")
RSFUNC(meddev, "realset-meddev")
RSFUNC(kurtosis, "realset-kurtosis")
RSFUNC(skewness, "realset-skewness")
RSFUNC(excess_kurtosis, "realset-excess-kurtosis")
RSFUNC(momental_skewness, "realset-momental-skewness")
RSFUNC(durbin_watson, "realset-durbin-watson")
RSFUNC(AR1_independence_z, "realset-AR1-independence-z")

GOOSE_PROC(autocorr, "realset-autocorr",
	   2,0,0, (SCM scm_rs, SCM scm_lag))
{
  SCM_ASSERT(REALSET_P(scm_rs), scm_rs, SCM_ARG1, str_autocorr);
  SCM_ASSERT(gh_exact_p(scm_lag), scm_lag, SCM_ARG2, str_autocorr);

  double x = 0;
  try {
    RealSet* rs = scm2realset(scm_rs);
    int lag = gh_scm2int(scm_lag);
    SCM_ASSERT(lag>=0 && lag<(int)rs->size()-3, scm_lag,
	       SCM_OUTOFRANGE, str_autocorr);
    x = goose_autocorr(*rs,lag);
  } GOOSE_CATCH(str_autocorr);

  return gh_double2scm(x);
}

GOOSE_PROC(autocorr_z, "realset-autocorr-z",
	   2,0,0, (SCM scm_rs, SCM scm_lag))
{
  SCM_ASSERT(REALSET_P(scm_rs), scm_rs, SCM_ARG1, str_autocorr_z);
  SCM_ASSERT(gh_exact_p(scm_lag), scm_lag, SCM_ARG2, str_autocorr_z);

  double x = 0;
  try {