plot.cpp

//**************************************************************************//
//* FRIDA: flexible rapid interactive data analysis                        *//
//* plot: plot on-line files                                               *//
//* (C) Joachim Wuttke 1990-, v2(C++) 2001-                                *//
//* http://apps.jcns.fz-juelich.de/frida                                   *//
//**************************************************************************//

#include <iostream>
#include <vector>
#include <string>

#include <read-plus/readln.hpp>

#include "defs.h"
#include "list.h"
#include "olf.h"
#include "mem.h"
#include "jsel.h"
#include "zentry.h"
#include "dualplot.h"
#include "plot.h"
#include "curve.h"

//**************************************************************************//
//*  Subroutines: determine data ranges                                    *//
//**************************************************************************//

//! Loop over selected data files to determine x range.

void determine_Xrange( CPlot* plot, CList* JSel )
{
    double inf, sup;
    NOlm::IterateD fiter;
    const vector<double> *vx;
    bool first = true;
    POld fd;
    while((fd=fiter())) {
        JSel->evaluate( 0, fd->nJ()-1 );
        for ( uint iv=0; iv<JSel->size(); ++iv ) {
            vx = &( fd->VS( JSel->V[iv] )->x );
            for ( uint i=0; i<vx->size(); ++i ){
                if ( plot->X.logflag && (*vx)[i]<=0 )
                    continue;
                if ( first ){
                    inf = (*vx)[i];
                    sup = (*vx)[i];
                    first = false;
                } else {
                    if ( (*vx)[i]<inf ) inf = (*vx)[i];
                    if ( (*vx)[i]>sup ) sup = (*vx)[i];
                }
            }
        }
    }
    if ( first )
        throw "x range is empty";
    plot->X.setRoundedLimits( "x", inf, sup );
}

//! Loop over selected files to determine y range for current x range.

void determine_Yrange( CPlot* plot, CList* JSel )
{
    double inf, sup;
    bool first=true;
    NOlm::IterateO iter;
    POlo f;
    while( (f=iter()) ){
        uint k = iter.k();
        POlc fc = P2C( f );
        if ( fc && fc->kconv!=-1 )
            continue;
        POld fd = P2D( f );
        JSel->evaluate( 0, f->nJ()-1 );
        for ( uint iv=0; iv<JSel->size(); ++iv ){ 
            uint j = JSel->V[iv];
            PSpec s;
            if ( fd ){
                s = fd->VS(j);
            } else {
                s = PSpec( new CSpec ); 
                plot->X.set_xgrid( s->x, 97 );
                fc->curve_val_vec( &(s->y), s->x, k, j );
            }
            for ( uint i=0; i<s->size(); i++ ) {
                if ( !plot->X.contains( s->x[i] ) )
                    continue;
                if ( plot->Y.logflag && s->y[i]<=0 )
                    continue;
                if (first) {
                    inf=s->y[i], sup=s->y[i];
                    first = false;
                } else {
                    if (s->y[i]<inf) inf=s->y[i];
                    if (s->y[i]>sup) sup=s->y[i];
                }
            }
        }
    }
    if ( first )
        throw "y range is empty";
    plot->Y.setRoundedLimits( "y", inf, sup );
}


//**************************************************************************//
//*  Subroutines: plot one scan                                            *//
//**************************************************************************//

//! Plot scan j of data file fd; return number of points plotted.

int plot_data( CPlot* plot, POld fd, uint k, uint j, int pstyle )
{
    const PSpec s = fd->VS(j);
    uint n=s->x.size();
    if ( n!=s->y.size() )
        throw "BUG: plot: x.size<>y.size";
    // Filter: x,y -> xp,yp if inside frame
    uint np=0, nxl=0, nxh=0, nyl=0, nyh=0;
    vector<double> xp, yp, dyp;
    for ( uint i=0; i<n; i++ ) {
        double x = s->x[i];
        if ( !plot->X.contains(x) ) {
            if ( x<=plot->X.inf ){
                x = plot->X.inf;
                nxl++;
            }
            if ( x>=plot->X.sup ){
                x = plot->X.sup;
                nxh++;
            }
            if ( !plot->X.force )
                continue;
        }
        double y = s->y[i];
        if ( !plot->Y.contains(y) ) {
            if ( y<=plot->Y.inf ){
                y = plot->Y.inf;
                nyl++;
            }
            if ( y>=plot->Y.sup ){
                y = plot->Y.sup;
                nyh++;
            }
            if ( !plot->Y.force )
                continue;
        }
        if ( np==plot->maxpoints && np<n ) {
            cout << "reached maxpoints at " << s->x[i] << "\n";
        }
        xp.push_back( x );
        yp.push_back( y );
        if ( s->dy.size() )
            dyp.push_back( s->dy[i] );
        np = xp.size();
        if ( np > plot->maxpoints ) // continue in large steps
            i += plot->maxpoints;
    }
    if ( np==0 ) {
        cout << "file " << k << " spec " << j << " all " << n <<
            " points outside plot range:\n";
        if ( nxl )
            cout << "  " << nxl << " points < xmin\n";
        if ( nxh )
            cout << "  " << nxh << " points > xmax\n";
        if ( nyl )
            cout << "  " << nyl << " points < ymin\n";
        if ( nyh )
            cout << "  " << nyh << " points > ymax\n";
    } else {
        plot->addSpec(
            false, true, pstyle, xp, yp, dyp, fd->V[j]->z,
            fd->xco.str(), fd->yco.str(),
            "data file " + S(k) + " spectrum "+S(j) );
    }
    return np;
}

//! Plot scan j of convolved curve file fc; return number of points plotted.

int plot_curve_convolved( CPlot* plot, POlc fc, uint k, uint j, int cstyle )
{
    vector<double> novec;
    uint kconv, jconv;
    PSpec sconv;
    NCurveFile::setConv( &sconv, &kconv, &jconv, k, j );
    vector<double> xp, yp;
    for ( uint i=0; i<sconv->size(); ++i ){
        double x = sconv->x[i];
        if ( plot->X.contains( x ) )
            xp.push_back( x );
    }
    fc->curve_val_vec( &yp, xp, k, j );
    vector<double> xo, yo;
    for ( uint i=0; i<xp.size(); ++i ){
        if ( !( plot->X.contains( xp[i] ) &&
                plot->Y.contains( yp[i] )   ) )
            continue;
        xo.push_back( xp[i] );
        yo.push_back( yp[i] );
    }
    if ( xo.size()==0 ){
        cout << "curve k="<<S(k)<<", j="<<S(j)<<
            " has no points in plot window\n";
    }
    plot->addSpec( true, true, cstyle, xo, yo, novec, fc->V[j]->z,
                   fc->xco.str(), fc->yco.str(),
                   "curve file "+S(k)+" spectrum "+S(j) );
    return xo.size();
}

//! Plot scan j of non-convolved curve file fc on equidistant grid;
//! return number of points plotted.

int plot_curve_equidist( CPlot* plot, POlc fc, uint k, uint j, int cstyle )
{
    vector<double> novec;
    // equidistant grid
    
    uint npts = plot->equipoints;
    vector<double> xc, yc, xp, yp;
    plot->X.set_xgrid( xc, npts );
    fc->curve_val_vec( &yc, xc, k, j );
    for ( uint i=0; i<npts; ++i ) {
        if ( plot->X.contains(xc[i]) && plot->Y.contains(yc[i]) ) {
            xp.push_back( xc[i] );
            yp.push_back( yc[i] );
        }
    }
    plot->addSpec( true, true, cstyle, xp, yp, novec, fc->V[j]->z,
                   fc->xco.str(), fc->yco.str(),
                   "curve file "+S(k)+" spectrum "+S(j) );
    return xp.size();
}

//! Plot scan j of non-convolved curve file fc, refining the grid where
//! necessary; return number of points plotted.

int plot_curve_refine( CPlot* plot, POlc fc, uint k, uint j, int cstyle )
{
    // refinement cares to:
    // - stop/start line when leaving/reentering/crossing(TODO) the frame
    // - TODO interpolate when crossing the horizontal border
    // - interpolate at discontinuities     
    
    vector<double> xc, yc, xn, yn;
    vector<double> novec;
                
    // start with equidistant grid:
    uint npts = plot->equipoints;
    plot->X.set_xgrid( xc, npts );

    // refinement loop:
    for ( int iref=0; ; ++iref ) {
        // cout << "DEB iref=" << iref << "\n";

        // evaluate curve for grid xc:
        fc->curve_val_vec( &yc, xc, k, j );
        // for ( uint i=0; i<xc.size(); ++i )
        // printf ("DEB xc yc %12.6g %12.6g\n", xc[i], yc[i] );
        // cout << "DEB evaluation -> " << xc.size() << "\n";

        // merge xc,yc and xn,yn:
        if ( iref==0 ) {
            xn = xc; yn = yc;
        } else {
            vector<double> xa, ya;
            uint ic=0, in=0;
            for ( ; ic<xc.size() && in<xn.size(); ) {
                if ( plot->X.close_enough( xc[ic], xn[ic], 1e-7 ) ){
                    ++in; // skip new point
                } else if ( xc[ic]<xn[in] ) {
                    xa.push_back( xc[ic] );
                    ya.push_back( yc[ic] );
                    ++ic;
                } else {
                    xa.push_back( xn[in] );
                    ya.push_back( yn[in] );
                    ++in;
                }
            }
            for ( ; ic<xc.size(); ++ic ) {
                xa.push_back( xc[ic] );
                ya.push_back( yc[ic] );
            }
            for ( ; in<xn.size(); ++in ) {
                xa.push_back( xn[in] );
                ya.push_back( yn[in] );
            }
            xn = xa; yn = ya;
        }
        // cout << "DEB merger -> " << xn.size() << "\n";
        // for ( uint i=0; i<xn.size(); ++i )
        // printf ("DEB xn yn %12.6g %12.6g\n", xn[i], yn[i] );

        // end of loop
        if ( iref>=10 || xn.size()>=plot->maxpoints )
            break;

        // eliminate points outside y range (except border points):
        vector<double> xa, ya;
        xa.push_back( xn[0] );
        ya.push_back( yn[0] );
        for ( uint i=1; i<xn.size()-1; ++i) {
        new_segment:
            xa.push_back( xn[i] );
            ya.push_back( yn[i] );
            if ( !plot->Y.contains(yn[i]) ) {
                for ( uint ii=i+1; ii<xn.size()-1; ++ii ) {
                    if ( plot->Y.contains(yn[ii]) ) {
                        if ( ii-1>i ) {
                            xa.push_back( xn[ii-1] );
                            ya.push_back( yn[ii-1] );
                        }
                        i = ii;
                        goto new_segment;
                    }
                }
            }
        }
        xa.push_back( xn.back() );
        ya.push_back( yn.back() );
        xn = xa; yn = ya;
        // cout << "DEB within frame -> " << xn.size() << "\n";
        // for ( uint i=0; i<xn.size(); ++i )
        // printf ("DEB xy %12.6g %12.6g\n", xn[i], yn[i] );

        // set additional base points:
        bool insert_next = false;
        xc.clear();
        for ( uint i=1; i<xn.size()-1; ++i ) {
            if ( i<xn.size()-2 ) {
                double yi = ( (xn[i+1]-xn[i])*yn[i-1] +
                              (xn[i]-xn[i-1])*yn[i+1] ) /
                    (xn[i+1]-xn[i-1]);
                if ( !plot->Y.close_enough( yn[i], yi, 0.005 ) ) {
                    xc.push_back( (xn[i-1] + xn[i])/2 );
                    // cout << "DEB disc " << xc.back() << "\n";
                    insert_next = true;
                    continue;
                }
            }
            if ( insert_next ) {
                xc.push_back( (xn[i-1] + xn[i])/2 );
                // cout << "DEB next " << xc.back() << "\n";
                insert_next = false;
                continue;
            }
            if ( ( plot->Y.contains(yc[i]) &&
                   !plot->Y.contains(yc[i-1]) ) ||
                 ( !plot->Y.contains(yc[i]) &&
                   plot->Y.contains(yc[i-1]) ) ) {
                // cout << "DEB edge " << xc.back() << "\n";
                xc.push_back( (xn[i-1] + xn[i])/2 );
                continue;
            }
        }
        // temporary check
        if( xc.size() )
            for ( uint i=0; i<xc.size()-1; ++i )
                if ( xc[i+1]<=xc[i] )
                    throw "BUG: new base points not sorted at " +
                        S(i) + ": " + S(xc[i]) + " " +
                        S(xc[i+1]);

    }

    // divide into segments and plot:
    vector<double> xa, ya;
    bool first_seg = true;
    npts = 0;
    for ( uint i=0; i<xn.size(); ++i ) {
        if ( plot->Y.contains(yn[i]) ) {
            xa.push_back( xn[i] );
            ya.push_back( yn[i] );
        } 
        if ( !plot->Y.contains(yn[i]) || i==xn.size()-1 ) {
            npts += xa.size();
            if ( xa.size() ){
                // cout << "DEB segment -> " << xa.size() << "\n";
                plot->addSpec(
                    true, first_seg, cstyle, xa, ya,
                    novec, fc->V[j]->z, fc->xco.str(), fc->yco.str(),
                    "curve file "+S(k)+ " spectrum "+S(j) );
                xa.clear();
                ya.clear();
                first_seg = false;
            }
        }
    }
    return npts;
}


//**************************************************************************//
//*  Plot routine is accessible from the main menu                         *//
//**************************************************************************//

//! Plot set of spectra, or add to open plot.

void NPlot::Plot( class CPlot *plot, bool add, const string& mode )
{
    NOlm::IterateO fiter;

    static CList JSel;
    static int pstyle=1, cstyle=1;
    static string jSel = "";

    if        ( mode=="ask" ) {
        JSelAsk( "Plot which spectra", &jSel, &JSel );
    } else if ( mode == "next" ) {
        JSel.increment( +1, 0, NOlm::nJ_max()-1 );
        cout << "plotting spectra " << JSel << "\n";
        jSel = JSel.str();
    } else if ( mode == "prev" ) {
        JSel.increment( -1, 0, NOlm::nJ_max()-1 );
        cout << "plotting spectra " << JSel << "\n";
        jSel = JSel.str();
    }

    if ( !add ) { // prepare frame (ranges, labels)
        pstyle = 1;
        cstyle = 1;
        if (!plot->X.finite())
            determine_Xrange( plot, &JSel );
        if (!plot->Y.finite())
            determine_Yrange( plot, &JSel );

        // build label: (zur Zeit nur vom ersten File; definiere +=)
        CCoord xCo = NOlm::sel_first()->xco;
        CCoord yCo = NOlm::sel_first()->yco;

        // draw new frame:
        plot->clearFrame();
        plot->plotFrame( xCo.ps_str(), yCo.ps_str() );

    }

    // plot:
    while( POlo f = fiter() ){
        uint k = fiter.k();
        POld fd = P2D( f );
        POlc fc = P2C( f );

        plot->docTxLine( f->name );
        for (uint i=0; i<f->lDoc.size(); i++)
            plot->docTxLine( "  " + f->lDoc[i] );
        if ( fc ) {
            vector<string> vs = fc->infoFile();
            for (uint i=0; i<vs.size(); i++)
                plot->docTxLine( "  " + vs[i] );
            plot->docTxLine( "  " + fc->infoScanHeader() );
        }

        JSel.evaluate( 0, f->nJ()-1 );
        for ( uint iv=0; iv<JSel.size(); ++iv){
            uint j = JSel.V[iv];

            if        ( fd ) {
                if ( !plot_data( plot, fd, k, j, pstyle ) )
                    continue;

            } else if ( fc ) {
                if        ( fc->kconv!=-1 ) {
                    if ( !plot_curve_convolved( plot, fc, k, j, cstyle ) )
                        continue;
                } else if ( !plot->refine ) {
                    if ( !plot_curve_equidist( plot, fc, k, j, cstyle ) )
                        continue;
                } else {
                    if ( !plot_curve_refine( plot, fc, k, j, cstyle ) )
                        continue;
                }

            } else {
                throw "PROGRAM ERROR plot: unexpected file type";
            }

            if        ( fd ) {
                plot->docPtTxLine( "  " + f->infoLine( j ), pstyle );
                ++pstyle;
            } else if ( fc ) {
                plot->docCvTxLine( "  " + f->infoLine( j ), cstyle );
                ++cstyle;
            }
        }
    }
    plot->showSpecs();
}