/****************************************************************************/
/* (C) Copyright 1993 Universidade Estadual de Campinas (UNICAMP)           */
/*                    Campinas, SP, Brazil                                  */
/*                                                                          */
/* This file can be freely distributed, modified, and used for any          */
/*   non-commercial purpose, provided that this copyright and authorship    */
/*   notice be included in any copy or derived version of this file.        */
/*                                                                          */
/* DISCLAIMER: This software is offered ``as is'', without any guarantee    */
/*   as to fitness for any particular purpose.  Neither the copyright       */
/*   holder nor the authors or their employers can be held responsible for  */
/*   any damages that may result from its use.                              */
/****************************************************************************/

/* See aaquad.h */

#include "aaquad.h"
#include "fltzeros2.h"
#include <aa.h>
#include <ia.h>
#include <flt.h>
#include <js.h>
#include <ps.h>
#include <math.h>
#include <ioprotos.h>
#include <stdio.h>
#include <malloc.h>

/*** PROTOTYPES FOR INTERNAL ROUTINES ***/

void aaquad_plot_aa_zeros(
    FILE *psfile,
    int epsformat,
    AAP ff (AAP x, AAP y),
    Interval xd,
    Interval yd,
    int n
  );
  /* 
    Plots roots of ff(xf,yf) = 0, using AA arithmetic.
    The input intervals are nxn cells in xd x yd. */

char *aaquad_format_parms(
    Interval xd,
    Interval yd,
    int n
  );
  /* Formats the arguments into a string */
  
void aaquad_report_stats(
    FILE *psfile,
    int epsformat, 
    char *arith,    /* "IA" or "AA" */
    int nevals, 
    int ncells
  );
  /* Prints evaluation report, and adds caption lines if not $epsformat$ */
  
void aaquad_aa_zeros_aux(
    FILE *psfile,
    AAP ff (AAP x, AAP y),
    Interval xd,
    Interval yd,
    int n,
    int *nevals,
    int *ncells,
    int ixlo, int mx,
    int iylo, int my
  );
  /* Called by aaquad_plot_aa_zeros */

void aaquad_plot_ia_zeros(
    FILE *psfile,
    int epsformat,
    Interval fv (Interval x, Interval y),
    Interval xd,
    Interval yd,
    int n
  );
  /* 
    Plots zeros of fv(xv,yv) = 0, using ordinary interval arithmetic.
    The input inervals are nxn cells in xd x yd. */

void aaquad_ia_zeros_aux(
    FILE *psfile,
    Interval fv (Interval x, Interval y),
    Interval xd,
    Interval yd,
    int n,
    int *nevals,
    int *ncells,
    int ixlo, int mx,
    int iylo, int my
  );
  /* Called by aaquad_plot_ia_zeros */

/*** IMPLEMENTATIONS ***/

void aaquad_plots(
    char *fileprefix,
    int epsformat,
    char *title,
    Float f(Float x, Float y),
    Interval fv(Interval x, Interval y),
    AAP ff(AAP x, AAP y),
    Interval xd,
    Interval yd,
    int n,
    int m
  )
  {
    char *parmstr = aaquad_format_parms(xd, yd, n);
    FILE *psfile;
    Float wx = xd.hi - xd.lo;
    Float wy = yd.hi - yd.lo;
    Float wm = FMAX(wx, wy);
    double scale = 6.0 * 72.0 / wm;
    double xc = 4.25 * 72.0;
    double hmin = xc - scale * wx / 2.0;
    double hmax = xc + scale * wx / 2.0;
    double yc = 6.50 *72.0;
    double vmin = yc - scale * wy / 2.0;
    double vmax = yc + scale * wy / 2.0;

    /*** Ordinary interval arithmetic plot ***/
    
    if (epsformat)
      {
        psfile = open_write(txtcat(fileprefix, "-ia.eps"));
        ps_begin_figure(
	  psfile,
	  xd.lo, xd.hi, yd.lo, yd.hi,
	  hmin,  hmax,  vmin,  vmax,
	  n, n
	);
      }
    else
      {
        psfile = open_write(txtcat(fileprefix, ".ps"));
        ps_begin_document(psfile);
        ps_begin_page(
	  psfile, 1,
	  xd.lo, xd.hi, yd.lo, yd.hi,
	  hmin,  hmax,  vmin,  vmax,
	  n, n
	);
        ps_add_caption(psfile, title);
        ps_add_caption(psfile, "");
        ps_add_caption(psfile, parmstr);
        ps_add_caption(psfile, "");
        ps_add_caption(psfile, "Ordinary interval arithmetic");
      }

    aaquad_plot_ia_zeros(psfile, epsformat, fv, xd, yd, n);
    fltzeros2_plot(psfile, f, xd, yd, m);
    ps_draw_frame(psfile);

    if (epsformat)
      { 
        ps_end_figure(psfile);
        fclose(psfile);
      }
    else
      {
        ps_end_page(psfile);
      }
      
    /*** Affine arithmetic plot ***/

    if (epsformat)
      {
        psfile = open_write(txtcat(fileprefix, "-aa.eps"));
	ps_begin_figure(
	  psfile,
	  xd.lo, xd.hi, yd.lo, yd.hi,
	  hmin,  hmax,  vmin,  vmax,
	  n, n
	);
      }
    else
      {
        ps_begin_page(
	  psfile, 2,
	  xd.lo, xd.hi, yd.lo, yd.hi,
	  hmin,  hmax,  vmin,  vmax,
	  n, n
	);
	ps_add_caption(psfile, title);
	ps_add_caption(psfile, "");
	ps_add_caption(psfile, parmstr);
	ps_add_caption(psfile, "");
	ps_add_caption(psfile, "Affine arithmetic");
      }

    aaquad_plot_aa_zeros(psfile, epsformat, ff, xd, yd, n);
    fltzeros2_plot(psfile, f, xd, yd, m);
    ps_draw_frame(psfile);

    if (epsformat) 
      {
        ps_end_figure(psfile);
        fclose(psfile);
      }
    else
      {
        ps_end_page(psfile);
        ps_end_document(psfile, 2);
        fclose(psfile);
      }
  }

char *aaquad_format_parms(
    Interval xd,
    Interval yd,
    int n
  )
  {
    char *s = (char *) malloc(100);
    sprintf(s, 
      "x in [%f _ %f]  y in [%f _ %f]  %d intervals",
      xd.lo, xd.hi, yd.lo, yd.hi, n
    );
    return(s);
  }

void aaquad_plot_aa_zeros(
    FILE *psfile,
    int epsformat,
    AAP ff (AAP x, AAP y),
    Interval xd,
    Interval yd,
    int n
  )
  { int nevals = 0;
    int ncells = 0;

    fprintf(stderr, "\nenter aaquad_plot_aa_zeros...\n");

    ps_begin_section(psfile, "Plot of AA-arithmetic zeros");

    aaquad_aa_zeros_aux(
      psfile, ff, xd, yd, n, &nevals, &ncells,
      0, n, 0, n
    );

    ps_end_section(psfile);

    fprintf(stderr, " done\n");
    
    aaquad_report_stats(psfile, epsformat, "AA", nevals, ncells);

    fprintf(stderr, "exit aaquad_plot_aa_zeros.\n");
  }

void aaquad_aa_zeros_aux(
    FILE *psfile,
    AAP ff (AAP x, AAP y),
    Interval xd,
    Interval yd,
    int n,
    int *nevals,
    int *ncells,
    int ixlo, int mx,
    int iylo, int my
  )
  {
    Interval xv, yv, zv;
    AAP xf, yf, zf;
    double gray = 0.75;

    MemP frame = aa_top();

    fprintf(stderr, "(");

    ROUND_DOWN;
    xv.lo = xd.lo + ((xd.hi - xd.lo)*ixlo)/n;
    yv.lo = yd.lo + ((yd.hi - yd.lo)*iylo)/n;

    ROUND_UP;
    xv.hi = xd.lo + ((xd.hi - xd.lo)*(ixlo+mx))/n;
    yv.hi = yd.lo + ((yd.hi - yd.lo)*(iylo+my))/n;

    xf = aa_from_interval(xv);
    yf = aa_from_interval(yv);

    zf = ff(xf, yf);
    zv = aa_range(zf);
    (*nevals)++;

    aa_flush(frame);

    ROUND_NEAR;
    if ((zv.lo > Zero) || (zv.hi < Zero))
      { ps_draw_rectangle(psfile, xv.lo, xv.hi, yv.lo, yv.hi); }
    else if ((mx == 1) && (my == 1))
      { (*ncells)++;
        ps_fill_and_draw_rectangle(psfile, xv.lo, xv.hi, yv.lo, yv.hi, gray);
      }
    else if (mx >= my)
      { int mx2 = mx / 2;
        aaquad_aa_zeros_aux(
          psfile, ff, xd, yd, n, nevals, ncells,
          ixlo, mx2, iylo, my
        );
        aaquad_aa_zeros_aux(
          psfile, ff, xd, yd, n, nevals, ncells,
          ixlo + mx2, mx - mx2, iylo, my
        );
      }
    else /* if (my > mx) */
      { int my2 = my / 2;
        aaquad_aa_zeros_aux(
          psfile, ff, xd, yd, n, nevals, ncells,
          ixlo, mx, iylo, my2
        );
        aaquad_aa_zeros_aux(
          psfile, ff, xd, yd, n, nevals, ncells,
          ixlo, mx, iylo + my2, my - my2
        );
      }

    fprintf(stderr, ")");
  }

void aaquad_plot_ia_zeros(
    FILE *psfile,
    int epsformat,
    Interval fv (Interval x, Interval y),
    Interval xd,
    Interval yd,
    int n
  )
  { int nevals = 0;
    int ncells = 0;

    fprintf(stderr, "\nenter aaquad_plot_ia_zeros...\n");

    ps_begin_section(psfile, "Plot of interval-arithmetic zeros");

    aaquad_ia_zeros_aux(
      psfile, fv, xd, yd, n, &nevals, &ncells,
      0, n, 0, n
    );

    ps_end_section(psfile);

    fprintf(stderr, "\n");
    
    aaquad_report_stats(psfile, epsformat, "IA", nevals, ncells);

    fprintf(stderr, "exit aaquad_plot_ia_zeros\n");
  }

void aaquad_report_stats(
    FILE *psfile, 
    int epsformat, 
    char *arith, 
    int nevals, 
    int ncells
  )
  {
    char sevals[80];

    fprintf(stderr, "%d %s function evaluations\n", nevals, arith);
    if (! epsformat)
      { sprintf(sevals, "%d %s function evaluations", nevals, arith);
        ps_add_caption(psfile, sevals);
      }

    fprintf(stderr, "%d cells retained\n", ncells);
    if (! epsformat)
      { sprintf(sevals, "%d cells retained", ncells);
        ps_add_caption(psfile, sevals);
      }
  }

void aaquad_ia_zeros_aux(
    FILE *psfile,
    Interval fv (Interval x, Interval y),
    Interval xd,
    Interval yd,
    int n,
    int *nevals,
    int *ncells,
    int ixlo, int mx,
    int iylo, int my
  )
  {
    Interval xv, yv, zv;
    double gray = 0.75;

    fprintf(stderr, "(");

    ROUND_DOWN;
    xv.lo = xd.lo + ((xd.hi - xd.lo)*ixlo)/n;
    yv.lo = yd.lo + ((yd.hi - yd.lo)*iylo)/n;

    ROUND_UP;
    xv.hi = xd.lo + ((xd.hi - xd.lo)*(ixlo+mx))/n;
    yv.hi = yd.lo + ((yd.hi - yd.lo)*(iylo+my))/n;

    zv = fv(xv, yv);
    (*nevals)++;

    ROUND_NEAR;
    if ((zv.lo > Zero) || (zv.hi < Zero))
      { ps_draw_rectangle(psfile, xv.lo, xv.hi, yv.lo, yv.hi); }
    else if ((mx == 1) && (my == 1))
      { ps_fill_and_draw_rectangle(psfile, xv.lo, xv.hi, yv.lo, yv.hi, gray);
        (*ncells)++;
      }
    else if (mx >= my)
      { int mx2 = mx / 2;
        aaquad_ia_zeros_aux(
          psfile, fv, xd, yd, n, nevals, ncells,
          ixlo, mx2, iylo, my
        );
        aaquad_ia_zeros_aux(
          psfile, fv, xd, yd, n, nevals, ncells,
          ixlo + mx2, mx - mx2, iylo, my
        );
      }
    else /* if (my > mx) */
      { int my2 = my / 2;
        aaquad_ia_zeros_aux(
          psfile, fv, xd, yd, n, nevals, ncells,
          ixlo, mx, iylo, my2
        );
        aaquad_ia_zeros_aux(
          psfile, fv, xd, yd, n, nevals, ncells,
          ixlo, mx, iylo + my2, my - my2
        );
      }

    fprintf(stderr, ")");
  }