dualplot.cpp

//**************************************************************************//
//* FRIDA: fast reliable interactive data analysis                         *//
//* dualplot.cpp: different mechanisms for screen and paper output         *//
//* (C) Joachim Wuttke 1990-, v2(C++) 2001-                                *//
//* http://frida.sourceforge.net                                           *//
//**************************************************************************//

#include <stdlib.h>
#include <string.h>
#include <iostream>
#include <math.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>
#include <fcntl.h>

#include "mystd.h"
#include "plotaux.h"
#include "coord.h"
#include "asi.h"
#include "gar.h"
#include "dualplot.h"

#define PSMAX 10
#define LINSIZ 80


//**************************************************************************//
//*   CRange                                                               *//
//**************************************************************************//

CRange::CRange(void)
{
	// initialize as infinite
	flag = 1;
	inf = 0;
	sup = 0;
	valid = 1;
}
	
CRange::CRange( double infin, double supin )
{
	flag = 0;
	inf = infin;
	sup = supin;
	test();
}	

CRange::CRange(int flagin)
{
	flag = flagin;
	inf = 0;
	sup = 0;
	valid = 1;
}
	

char RangeDecode( const string& resp, CRange *rval );

CRange::CRange( const string& text )
{
	CRange rval;
	char ok;
	ok = RangeDecode( text, &rval );
	if (ok!='v') {
		flag = 1;
		inf = 0;
		sup = 0;
		valid = 0;
	} else {
		*this = rval;
	}
}

void CRange::Round(int logflag, double relmargin, double digits)
{
	if (!test())	
		throw "Range::Round: illegal params: " + str();

	if        (inf==sup) {
		*this = CRange(inf-1., sup+1.);
	} else if (finite()) {
		double margin = relmargin * (sup - inf);
		double mydigits = digits;
		double inf2, sup2;
		do {
			if (mydigits>=8.5) return; // do not round
			sup2 = mystd::round_decimal( sup+margin, mydigits );
			inf2 = mystd::round_decimal( inf-margin, mydigits );
			if(sup<0 && sup2>0) sup2 = 0;
			if(inf2<0 && inf>0) inf2 = 0;
			mydigits += 0.5;
		} while (!((inf2<inf) && (sup<sup2)));
		*this = CRange(inf2, sup2);
	}
}	

bool CRange::infinite() const
{
	return (flag&1);
}

bool CRange::empty() const
{
	return (flag&2 || inf>sup || ((flag&4 || flag&8) && inf==sup));
}

bool CRange::finite() const
{
	return (!this->infinite() && !this->empty());
}

bool CRange::contained(double val)  const
{
	if (this->infinite())
		return 1;
	else if (this->empty())
		return 0;
	else
		return ((inf<val || (!(flag&4) && inf==val)) &&
			(sup>val || (!(flag&8) && sup==val)));
}

double CRange::relval(int logflag, double v) const
{
	if (flag&1 || flag&2) return 0;
	if (logflag) {
		if(inf<0 || v<0) return 0;
		return (double) 
		    (log(v/(double)inf) / log((double)sup/(double)inf));
	} else {
		return (v-inf) / (sup-inf);
	}
}

CRange CRange::operator*(const CRange& R)
{
	CRange Z;
	if        (infinite() || R.infinite()) {
		Z.flag = 1;
	} else if (R.empty()) {
		Z = *this;
	} else if (empty()) {
		Z = R;
	} else {
		// cout << " input: " << *this << " * " << R;
		Z = CRange(min(inf, R.inf), max(sup, R.sup));
		// Z.flag4 = ..
		// Z.flag8 = ..
		// cout << " output : " << Z << endl;
	}
	return Z;
}

CRange CRange::operator*=(const CRange& R)
{
	CRange Z = *this * R;
	return *this = Z;
}

bool CRange::test(void) 
{
	return valid = (inf<=sup || flag&1 || flag&2);
}

const char * CRange::helptext =
        "a floating-point range:\n"
        "   -> <x> <y> for a closed interval (x y)\n"
	"   -> -       for an empty range\n"
	"   -> *       for an infinite range\n";

char RangeDecode( const string& resp, CRange *rval )
{
	double r1, r2;
	if        ( resp=="" ) {
		return 'd';
	} else if ( resp=="*" ) {
		*rval = CRange(1);
		return 'v';
	} else if ( resp=="-" ) {
		*rval = CRange(2);
		return 'v';
	}
	string s1, s2;
        mystd::string_extract_word(string(resp), &s1, &s2);
	if (!mystd::any2dbl(s1,&r1) && !mystd::any2dbl(s2,&r2)) {
		*rval = CRange(r1, r2);
		if (!rval->test()) return 'h';
		return 'v';
	}
	return 'h';
}

CRange CRange::Ask( const string& quest )
{
	return Ask(quest, 0, CRange());
}

CRange CRange::Ask( const string& quest, const CRange& Def )
{
	return Ask(quest, 1, Def);
}

CRange CRange::Ask( const string& quest, const bool defflag,
	const CRange& Def)
{
	*this = myask( quest, defflag, Def, helptext, RangeDecode );
	return *this;
}

string CRange::str(void) const
{
	if      (!valid)
		return "<invalid range>";
	else if (finite()) 
		return ((flag&4) ? ")" : "(")
			+ strg(inf) + " " + strg(sup)
			+ ((flag&8) ? "(" : ")");
	else if (infinite())
		return "*";
	else
		return "-";
}

ostream& operator<< (ostream &s, const CRange& R)
{
	if      (!R.valid)
		return s << "<invalid range>";
	else if (R.finite()) 
		return s << ((R.flag&4) ? ")" : "(")
			 << R.inf << " " << R.sup 
			 << ((R.flag&8) ? "(" : ")");
	else if (R.infinite())
		return s << "*";
	else
		return s << "-";
}


//**************************************************************************//
//*  CAxis                                                                 *//
//**************************************************************************//

void CAxis::Ask(const string& quest)
{
    string def, in, in1, in2;
    CRange newR;
    while(1) {
        def = (log ? "log " : "") + R.str();
        // out = quest + " [" + def + "] ?"; so gehts leider(?) nicht 
        in = sask(quest, def);

        if        (in==def)
            return;

        mystd::string_extract_word(in, &in1, &in2);
        if        (in1=="l") {
            SetLog(!log);
            in = in2;
        } else if (in1=="i") {
            SetLog(0);
            in = in2;
        } else if (in1=="g") {
            SetLog(1);
            in = in2;
        } 

        if (in!="") {
            newR = CRange(in);
            if (!newR.valid)
                continue;
            if (log && newR.inf<=0) {
                printf("! log scale requires range above 0\n");
                continue;
            }
            if (!log && ( newR.inf<-3e38 || newR.sup>3e38 ) ) {
                printf( "! lin scale exceeds typical PS implementation range "
                        "(abs<3e38)\n" );
                continue;
            }
            R = newR;
            return;
        }
    }
}

void CAxis::SetLog(const uint val)
{
    if (log!=val) {
        log = val;
        CRange RX = AlternateR;
        AlternateR = R;
        R  = RX;
    }
}

double CAxis::pc(double v) // value -> plot-coordinate
{
    return PSMAX * R.relval(log, v);
}


//**************************************************************************//
//*  NPlot                                                                 *//
//**************************************************************************//

namespace NPlot {
    CAxis X, Y;

    uint maxpoints = 12000;

    int gp_fifo;
    void gp_write(string);
    int gp_fno;
    string gp_fnames;

    uint  ps_fnum; // file
    uint  ps_snum; // scan
    uint  ps_pnum; // scan with pstyle
    uint  ps_cnum; // scan with cstyle
    void ps_ticktack(uint ntack, double *tack, int ntpt, double *ticklim,
                     CAxis *A);
    vector<string> ps_Doc;
    vector<string> ps_accu; // future output is accumulated here
    char outlin[ LINSIZ ];
};

void NPlot::Open(void)
{
    // Start gnuplot. Use input redirection so that gnuplot receives
    // commands from a gp_fifo which is created here.

    string fn_fifo = string("/tmp/gnuplot-") + getenv( "LOGNAME" );

    system( ("rm -f "+fn_fifo).c_str() );
    if (mkfifo( fn_fifo.c_str(), 0666 )) {
        printf("SEVERE SYSTEM ERROR cannot make fifo %s"
               " - will not be able to print\n", fn_fifo.c_str() );
        return;
    }

    system( ("gnuplot -noraise < " + fn_fifo + " &").c_str() );

    // we use open instead of fopen or ofstream,
    // because we need non-blocking mode.
    if (!(gp_fifo = open(fn_fifo.c_str(), O_WRONLY))) {
        printf("SEVERE SYSTEM ERROR cannot open gp_fifo"
               " - will not be able to print\n");
        return;
    }
    fcntl(gp_fifo,F_SETFL,O_NONBLOCK);

    // gp_write(string("plot 1 2 3")); // test
    gp_write(string("set terminal x11"));

    ps_fnum=0;
}
		
void NPlot::Close(void)
{
}

void NPlot::Clear(void)
{
    // cout << "DEBUG dp Clear beg\n";
    gp_fno = 0;
    gp_fnames = string("");

    // new wups file:
    string cmd;
    ps_accu.clear();
    ps_accu.push_back( "\n%% output created by frida2\n");
    ps_snum = 0;
    ps_pnum = 0;
    ps_cnum = 0;
    ps_Doc.clear();
}

void NPlot::Multiask(void)
{
    maxpoints = iask("Max # points in plot", maxpoints);
    X.force = bask("Force x points into frame", X.force);
    Y.force = bask("Force y points into frame", Y.force);
}

void NPlot::Box()
{
    gp_write( "set nologscale" );
    string whichlog="";
    if (X.log) whichlog += "x";
    if (Y.log) whichlog += "y";
    if (whichlog!="")
        gp_write( "set logscale " + whichlog );

    // wups:
    snprintf(outlin, LINSIZ, "\n%d %g %g xSetCoord\n", 
             X.log, X.R.inf, X.R.sup);
    ps_accu.push_back( outlin );
    snprintf(outlin, LINSIZ, "%d %g %g ySetCoord\n", 
             Y.log, Y.R.inf, Y.R.sup);
    ps_accu.push_back( outlin );
    snprintf(outlin, LINSIZ, "%% %d %g %g %d zSetCoord\n\n", 0, 0., 0., 0);
    ps_accu.push_back( outlin );

    int ntack, ntpt;
    double *tack, ticklim[2];

    ps_accu.push_back( "\n/xPlotFrame {" );
    plotaux::calc_ticks( X.log, X.R.inf, X.R.sup, 
                         &ntack, &tack, &ntpt, ticklim );
    ps_ticktack(ntack, tack, ntpt, ticklim, &X);
    free(tack);
    snprintf( outlin, LINSIZ-4, "   {(%s", X.C.ps_str().c_str() );
    strncat( outlin, ")}\n", LINSIZ );
    ps_accu.push_back( outlin );
    ps_accu.push_back( "      0 10   0  0     0  90 "
                       "OneAxx Axx Tic xTacL xNumL %% low x axis\n" );
    ps_accu.push_back( "      0 10   0 10     0 270 "
                       "OneAxx Axx Tic xTacH       %% top x axis\n" );
    ps_accu.push_back( "      xCL\n" );
    ps_accu.push_back( "} def\n" );
    
    ps_accu.push_back( "\n/yPlotFrame {" );
    plotaux::calc_ticks( Y.log, Y.R.inf, Y.R.sup, 
                         &ntack, &tack, &ntpt, ticklim );
    ps_ticktack(ntack, tack, ntpt, ticklim, &Y);
    free(tack);
    snprintf( outlin, LINSIZ, "   {(%s", Y.C.ps_str().c_str() );
    strncat( outlin, ")}\n", LINSIZ );
    ps_accu.push_back( outlin );
    ps_accu.push_back( "      0 10   0  0    90   0 "
                       "OneAxx Axx Tic yTacL yNumL %% left y axis\n" );
    ps_accu.push_back( "      0 10  10  0    90 180 "
                       "OneAxx Axx Tic yTacH % yNumH %% right yaxis\n" );
    ps_accu.push_back( "      yCL\n" );
    ps_accu.push_back( "} def\n" );
    ps_accu.push_back( "\n%% modeDD\nplotbefore\n" );
}

void NPlot::Line(const int lstyle, 
		 const vector<double> x, const vector<double> y,
		 const vector<double>* z,
		 const string xco, const string yco)
{
    double *xp, *yp;
    uint np=0, nxl=0, nxh=0, nyl=0, nyh=0;
    uint n=x.size();
    if (n!=y.size()) {
        printf("PROG ERR NPLot::Line x.size<>y.size\n");
        return;
    }
    xp = (double *) malloc(sizeof(double)*n);
    yp = (double *) malloc(sizeof(double)*n);
    for (uint i=0; i<n; i++) {
        if(!X.force && !X.R.contained(x[i])) {
            if (x[i]<=X.R.inf) nxl++;
            if (x[i]>=X.R.sup) nxh++;
            continue;
        }
        if(!Y.force && !Y.R.contained(y[i])) {
            if (y[i]<=Y.R.inf) nyl++;
            if (y[i]>=Y.R.sup) nyh++;
            continue;
        }
        if(np==maxpoints && np<n) {
            printf("reached maxpoints at %g\n", x[i]);
        }
        xp[np] = x[i];
        yp[np] = y[i];
        np++;
        if(np>maxpoints) i += maxpoints;
    }
    if (np==0) {
        printf("no points in data range:\n"
               "of %d points are\n %d left and %d right of the x range,"
               "\n %d below and %d above of the y range\n", 
               n, nxl, nxh, nyl, nyh);
        return;
    }

    // save to file:
    char gp_fnam[40];
    sprintf( gp_fnam, "/tmp/%s-%03d.gnu", getenv("LOGNAME"), gp_fno++);
    if (gp_fnames!="") gp_fnames += ", ";
    gp_fnames += string("\"") + gp_fnam + "\"" + 
//		string(" title \"") + xco + string (" -> ") 
//		                    + yco + string ("\"");
        " notitle";
    if (lstyle<0)
        gp_fnames += " with lines";
    FILE *gp_fd;
    if (!(gp_fd = fopen(gp_fnam, "w"))) {
        printf( "cannot save gnuplot data to %s", gp_fnam );
        return;
    } 
    for (uint i=0; i<np; i++)
        fprintf(gp_fd, "%20.12g %20.12g\n", xp[i], yp[i]);
    fclose(gp_fd);
	
    // plot command:
    char aux[80];
    sprintf( aux, "[%12.8g:%12.8g] [%12.8g:%12.8g] ", 
             X.R.inf, X.R.sup, Y.R.inf, Y.R.sup);
    gp_write(string("plot ") + aux + gp_fnames);

    snprintf( outlin, LINSIZ, "\n%3u [", ++ps_snum);
    ps_accu.push_back( outlin );
    for (uint i=0; i<z->size(); i++){
        snprintf( outlin, LINSIZ, " %12g", (*z)[i]);
        ps_accu.push_back( outlin );
    }
    snprintf( outlin, LINSIZ, " ] zValues\n");
    ps_accu.push_back( outlin );

    if (lstyle>=0)
        snprintf( outlin, LINSIZ, "%2d pstyle", ++ps_pnum);
    else
        snprintf( outlin, LINSIZ, "%2d cstyle", 1);
    ps_accu.push_back( outlin );
    snprintf( outlin, LINSIZ-2, " %% (%s -> %s)", xco.c_str(), yco.c_str());
    strncat( outlin, "\n", LINSIZ );
    ps_accu.push_back( outlin );
    for (uint i=0; i<np; i++) {
        snprintf( outlin, LINSIZ,
                  "%7.3f%7.3f%7.3f t%c %% %14.7g wx %14.7g wy\n",
                  X.pc(xp[i]), Y.pc(yp[i]), 0.0,
                  i==0 ? 'i' : i==np-1 ? 'f' : ' ', xp[i], yp[i] );
        ps_accu.push_back( outlin );
    }

    free(xp);
    free(yp);
}

void NPlot::Doc (vector<string> lDoc)
{
    for (uint i=0; i<lDoc.size(); i++)
        ps_Doc.push_back(lDoc[i]);
}

void NPlot::Save( bool full_outfile )
{
    string ps_outdir, ps_head, ps_dict;

    // read configuration parameters:
    ps_outdir = NRead::wofmac( "\\psdir" );
    ps_head = NRead::wofmac( "\\pshead" );
        if ( full_outfile )
            ps_dict = NRead::wofmac( "\\psdict" );

    // construct output file name:
    FILE *pssav;
    char outf[20];
    string flong, cmd;
    while(1) {
        if (ps_fnum>=999)
            throw string( "graph file number overflow" );
        sprintf(outf, "%sl%d.%s", ps_outdir.c_str(), ++ps_fnum, 
                full_outfile ? "ps" : "psa" );
        if (!(pssav = mystd::glob_fopen(outf, "", "", "r")))
            break; // legal exit
        fclose(pssav);
    }
    printf("save plot in %s\n", outf);

    // copy headers to output file:
    cmd = string("cat ") + ( full_outfile ? ps_dict : "" ) + " " +
        ps_head + " > " + outf + "\n";
    if (system(cmd.c_str())) {
        printf ("cannot copy to %s\n", outf);
        return;
    }
    
    // a redundant check of existence of the outfile:
    if (!(pssav = mystd::glob_fopen(outf, "", "", "r", &flong))) {
        printf("have not created file %s\n", outf);
        return;
    }
    fclose(pssav);

    // append specific output to output file:
    if (!(pssav = fopen(flong.c_str(), "a+"))) {
        printf ("cannot write (append) to  file %s\n", flong.c_str());
        return;
    }
    for( uint i=0; i<ps_accu.size(); ++i ){
        // fprintf does not work here because output line may contain "%"
        fwrite( ps_accu[i].c_str(), 1, ps_accu[i].size(), pssav );
    }

    // additional output (do not append this to ps_accu to allow
    // further incrementation of ps_accu):
    fprintf(pssav, "\nblack 0 -4 13 newlist\n");
    for (uint i=0; i<ps_Doc.size(); i++) 
        fprintf(pssav, "{(%s)} infline\n", ps_Doc[i].c_str());
	
    fprintf(pssav, 
            "\n{(%s)}  /filename exch def 10 -2.8 18 showfilename\n\n"
            "EndFrame\n", outf);

    // output completed:
    fclose(pssav);
}

void NPlot::Dialog(void)
{
    string cmd;
    cout << "entering mygnuplot - to leave, type q\n";
    while(1) {
        cmd = sask("mygnuplot> ");
        if (cmd.substr(0,1)=="q") return;
        gp_write(cmd);
    }
}

void NPlot::gp_write(string in)
{
    // cout << "DEBUG gp_write " << in << "\n";
    string out = in + "\n";
    write(gp_fifo, out.c_str(), out.size());
}

void NPlot::ps_ticktack(uint ntack, double *tack, int ntpt, double *ticklim,
			CAxis *A)
{
    uint i;
    ps_accu.push_back( "[\n" );
    if (A->log && ( tack[0]<1e-3 || tack[ntack-1]>1e3 )) {
        for (i=0; i<ntack; i++) {
            snprintf( outlin, LINSIZ, "   %9.6f {(10)(%d)sp()} %%{(%g)}\n", 
                      A->pc(tack[i]), lrint(log10(tack[i])),
                      (float) tack[i]);
            ps_accu.push_back( outlin );
        }
    } else {
        for (i=0; i<ntack; i++) {
            snprintf( outlin, LINSIZ, "   %9.6f {(%g)}\n", 
                    A->pc(tack[i]), (float) tack[i]);
            ps_accu.push_back( outlin );
        }
    }
    ps_accu.push_back( "   ] SetTacVec\n" );
    snprintf( outlin, LINSIZ, "   %g %g %d %d SetTicVec%s\n", 
            A->pc(ticklim[0]), A->pc(ticklim[1]), ntack+2, ntpt,
            (A->log? "Log" : "Lin"));
    ps_accu.push_back( outlin );
}