text:
/*
 * This file is part of the GROMACS molecular simulation package.
 *
 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
 * Copyright (c) 2001-2004, The GROMACS development team,
 * check out http://www.gromacs.org for more information.
 * Copyright (c) 2012,2013, by the GROMACS development team, led by
 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
 * others, as listed in the AUTHORS file in the top-level source
 * directory and at http://www.gromacs.org.
 *
 * GROMACS is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public License
 * as published by the Free Software Foundation; either version 2.1
 * of the License, or (at your option) any later version.
 *
 * GROMACS 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
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with GROMACS; if not, see
 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
 *
 * If you want to redistribute modifications to GROMACS, please
 * consider that scientific software is very special. Version
 * control is crucial - bugs must be traceable. We will be happy to
 * consider code for inclusion in the official distribution, but
 * derived work must not be called official GROMACS. Details are found
 * in the README & COPYING files - if they are missing, get the
 * official version at http://www.gromacs.org.
 *
 * To help us fund GROMACS development, we humbly ask that you cite
 * the research papers on the package. Check out http://www.gromacs.org.
 */
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include "statutil.h"
#include "sysstuff.h"
#include "typedefs.h"
#include "smalloc.h"
#include "macros.h"
#include "gmx_fatal.h"
#include "vec.h"
#include "copyrite.h"
#include "futil.h"
#include "readinp.h"
#include "statutil.h"
#include "txtdump.h"
#include "gstat.h"
#include "xvgr.h"
#include "physics.h"
#include "gmx_ana.h"

enum {
    epAuf, epEuf, epAfu, epEfu, epNR
};
enum {
    eqAif, eqEif, eqAfi, eqEfi, eqAui, eqEui, eqAiu, eqEiu, eqNR
};
static char *eep[epNR] = { "Af", "Ef", "Au", "Eu" };
static char *eeq[eqNR] = { "Aif", "Eif", "Afi", "Efi", "Aui", "Eui", "Aiu", "Eiu" };

typedef struct {
    int       nreplica;  /* Number of replicas in the calculation                   */
    int       nframe;    /* Number of time frames                                   */
    int       nstate;    /* Number of states the system can be in, e.g. F,I,U       */
    int       nparams;   /* Is 2, 4 or 8                                            */
    gmx_bool *bMask;     /* Determine whether this replica is part of the d2 comp.  */
    gmx_bool  bSum;
    gmx_bool  bDiscrete; /* Use either discrete folding (0/1) or a continuous       */
    /* criterion */
    int       nmask;     /* Number of replicas taken into account                   */
    real      dt;        /* Timestep between frames                                 */
    int       j0, j1;    /* Range of frames used in calculating delta               */
    real    **temp, **data, **data2;
    int     **state;     /* State index running from 0 (F) to nstate-1 (U)          */
    real    **beta, **fcalt, **icalt;
    real     *time, *sumft, *sumit, *sumfct, *sumict;
    real     *params;
    real     *d2_replica;
} t_remd_data;

#ifdef HAVE_LIBGSL
#include <gsl/gsl_multimin.h>

static char *itoa(int i)
{
    static char ptr[12];

    sprintf(ptr, "%d", i);
    return ptr;
}

static char *epnm(int nparams, int index)
{
    static char buf[32], from[8], to[8];
    int         nn, ni, ii;

    range_check(index, 0, nparams);
    if ((nparams == 2) || (nparams == 4))
    {
        return eep[index];
    }
    else if ((nparams > 4) && (nparams % 4 == 0))
    {
        return eeq[index];
    }
    else
    {
        gmx_fatal(FARGS, "Don't know how to handle %d parameters", nparams);
    }

    return NULL;
}

static gmx_bool bBack(t_remd_data *d)
{
    return (d->nparams > 2);
}

static real is_folded(t_remd_data *d, int irep, int jframe)
{
    if (d->state[irep][jframe] == 0)
    {
        return 1.0;
    }
    else
    {
        return 0.0;
    }
}

static real is_unfolded(t_remd_data *d, int irep, int jframe)
{
    if (d->state[irep][jframe] == d->nstate-1)
    {
        return 1.0;
    }
    else
    {
        return 0.0;
    }
}

static real is_intermediate(t_remd_data *d, int irep, int jframe)
{
    if ((d->state[irep][jframe] == 1) && (d->nstate > 2))
    {
        return 1.0;
    }
    else
    {
        return 0.0;
    }
}

static void integrate_dfdt(t_remd_data *d)
{
    int    i, j;
    double beta, ddf, ddi, df, db, fac, sumf, sumi, area;

    d->sumfct[0] = 0;
    d->sumict[0] = 0;
    for (i = 0; (i < d->nreplica); i++)
    {
        if (d->bMask[i])
        {
            if (d->bDiscrete)
            {
                ddf = 0.5*d->dt*is_folded(d, i, 0);
                ddi = 0.5*d->dt*is_intermediate(d, i, 0);
            }
            else
            {
                ddf = 0.5*d->dt*d->state[i][0];
                ddi = 0.0;
            }
            d->fcalt[i][0] = ddf;
            d->icalt[i][0] = ddi;
            d->sumfct[0]  += ddf;
            d->sumict[0]  += ddi;
        }
    }
    for (j = 1; (j < d->nframe); j++)
    {
        if (j == d->nframe-1)
        {
            fac = 0.5*d->dt;
        }
        else
        {
            fac = d->dt;
        }
        sumf = sumi = 0;
        for (i = 0; (i < d->nreplica); i++)
        {
            if (d->bMask[i])
            {
                beta = d->beta[i][j];
                if ((d->nstate <= 2) || d->bDiscrete)
                {
                    if (d->bDiscrete)
                    {
                        df = (d->params[epAuf]*exp(-beta*d->params[epEuf])*
                              is_unfolded(d, i, j));
                    }
                    else
                    {
                        area = (d->data2 ? d->data2[i][j] : 1.0);
                        df   =  area*d->params[epAuf]*exp(-beta*d->params[epEuf]);
                    }
                    if (bBack(d))
                    {
                        db = 0;
                        if (d->bDiscrete)
                        {
                            db = (d->params[epAfu]*exp(-beta*d->params[epEfu])*
                                  is_folded(d, i, j));
                        }
                        else
                        {
                            gmx_fatal(FARGS, "Back reaction not implemented with continuous");
                        }
                        ddf = fac*(df-db);
                    }
                    else
                    {
                        ddf = fac*df;
                    }
                    d->fcalt[i][j] = d->fcalt[i][j-1] + ddf;
                    sumf          += ddf;
                }
                else
                {
                    ddf = fac*((d->params[eqAif]*exp(-beta*d->params[eqEif])*
                                is_intermediate(d, i, j)) -
                               (d->params[eqAfi]*exp(-beta*d->params[eqEfi])*
                                is_folded(d, i, j)));
                    ddi = fac*((d->params[eqAui]*exp(-beta*d->params[eqEui])*
                                is_unfolded(d, i, j)) -
                               (d->params[eqAiu]*exp(-beta*d->params[eqEiu])*
                                is_intermediate(d, i, j)));
                    d->fcalt[i][j] = d->fcalt[i][j-1] + ddf;
                    d->icalt[i][j] = d->icalt[i][j-1] + ddi;
                    sumf          += ddf;
                    sumi          += ddi;
                }
            }
        }
        d->sumfct[j] = d->sumfct[j-1] + sumf;
        d->sumict[j] = d->sumict[j-1] + sumi;
    }
    if (debug)
    {
        fprintf(debug, "@type xy\n");
        for (j = 0; (j < d->nframe); j++)
        {
            fprintf(debug, "%8.3f  %12.5e\n", d->time[j], d->sumfct[j]);
        }
        fprintf(debug, "&\n");
    }
}

static void sum_ft(t_remd_data *d)
{
    int    i, j;
    double fac;

    for (j = 0; (j < d->nframe); j++)
    {
        d->sumft[j] = 0;
        d->sumit[j] = 0;
        if ((j == 0) || (j == d->nframe-1))
        {
            fac = d->dt*0.5;
        }
        else
        {
            fac = d->dt;
        }
        for (i = 0; (i < d->nreplica); i++)
        {
            if (d->bMask[i])
            {
                if (d->bDiscrete)
                {
                    d->sumft[j] += fac*is_folded(d, i, j);
                    d->sumit[j] += fac*is_intermediate(d, i, j);
                }
                else
                {
                    d->sumft[j] += fac*d->state[i][j];
                }
            }
        }
    }
}

static double calc_d2(t_remd_data *d)
{
    int    i, j;
    double dd2, d2 = 0, dr2, tmp;

    integrate_dfdt(d);

    if (d->bSum)
    {
        for (j = d->j0; (j < d->j1); j++)
        {
            if (d->bDiscrete)
            {
                d2  += sqr(d->sumft[j]-d->sumfct[j]);
                if (d->nstate > 2)
                {
                    d2 += sqr(d->sumit[j]-d->sumict[j]);
                }
            }
            else
            {
                d2  += sqr(d->sumft[j]-d->sumfct[j]);
            }
        }
    }
    else
    {
        for (i = 0; (i < d->nreplica); i++)
        {
            dr2 = 0;
            if (d->bMask[i])
            {
                for (j = d->j0; (j < d->j1); j++)
                {
                    tmp  = sqr(is_folded(d, i, j)-d->fcalt[i][j]);
                    d2  += tmp;
                    dr2 += tmp;
                    if (d->nstate > 2)
                    {
                        tmp  = sqr(is_intermediate(d, i, j)-d->icalt[i][j]);
                        d2  += tmp;
                        dr2 += tmp;
                    }
                }
                d->d2_replica[i] = dr2/(d->j1-d->j0);
            }
        }
    }
    dd2 = (d2/(d->j1-d->j0))/(d->bDiscrete ? d->nmask : 1);

    return dd2;
}

static double my_f(const gsl_vector *v, void *params)
{
    t_remd_data *d       = (t_remd_data *) params;
    double       penalty = 0;
    int          i;

    for (i = 0; (i < d->nparams); i++)
    {
        d->params[i] = gsl_vector_get(v, i);
        if (d->params[i] < 0)
        {
            penalty += 10;
        }
    }
    if (penalty > 0)
    {
        return penalty;
    }
    else
    {
        return calc_d2(d);
    }
}

static void optimize_remd_parameters(FILE *fp, t_remd_data *d, int maxiter,
                                     real tol)
{
    real   size, d2;
    int    iter   = 0;
    int    status = 0;
    int    i;

    const gsl_multimin_fminimizer_type *T;
    gsl_multimin_fminimizer            *s;

    gsl_vector                         *x, *dx;
    gsl_multimin_function               my_func;

    my_func.f      = &my_f;
    my_func.n      = d->nparams;
    my_func.params = (void *) d;

    /* Starting point */
    x = gsl_vector_alloc (my_func.n);
    for (i = 0; (i < my_func.n); i++)
    {
        gsl_vector_set (x, i, d->params[i]);
    }

    /* Step size, different for each of the parameters */
    dx = gsl_vector_alloc (my_func.n);
    for (i = 0; (i < my_func.n); i++)
    {
        gsl_vector_set (dx, i, 0.1*d->params[i]);
    }

    T = gsl_multimin_fminimizer_nmsimplex;
    s = gsl_multimin_fminimizer_alloc (T, my_func.n);

    gsl_multimin_fminimizer_set (s, &my_func, x, dx);
    gsl_vector_free (x);
    gsl_vector_free (dx);

    printf ("%5s", "Iter");
    for (i = 0; (i < my_func.n); i++)
    {
        printf(" %12s", epnm(my_func.n, i));
    }
    printf (" %12s %12s\n", "NM Size", "Chi2");

    do
    {
        iter++;
        status = gsl_multimin_fminimizer_iterate (s);

        if (status != 0)
        {
            gmx_fatal(FARGS, "Something went wrong in the iteration in minimizer %s",
                      gsl_multimin_fminimizer_name(s));
        }

        d2     = gsl_multimin_fminimizer_minimum(s);
        size   = gsl_multimin_fminimizer_size(s);
        status = gsl_multimin_test_size(size, tol);

        if (status == GSL_SUCCESS)
        {
            printf ("Minimum found using %s at:\n",
                    gsl_multimin_fminimizer_name(s));
        }

        printf ("%5d", iter);
        for (i = 0; (i < my_func.n); i++)
        {
            printf(" %12.4e", gsl_vector_get (s->x, i));
        }
        printf (" %12.4e %12.4e\n", size, d2);
    }
    while ((status == GSL_CONTINUE) && (iter < maxiter));

    gsl_multimin_fminimizer_free (s);
}

static void preprocess_remd(FILE *fp, t_remd_data *d, real cutoff, real tref,
                            real ucut, gmx_bool bBack, real Euf, real Efu,
                            real Ei, real t0, real t1, gmx_bool bSum, gmx_bool bDiscrete,
                            int nmult)
{
    int  i, j, ninter;
    real dd, tau_f, tau_u;

    ninter = (ucut > cutoff) ? 1 : 0;
    if (ninter && (ucut <= cutoff))
    {
        gmx_fatal(FARGS, "You have requested an intermediate but the cutoff for intermediates %f is smaller than the normal cutoff(%f)", ucut, cutoff);
    }

    if (!bBack)
    {
        d->nparams = 2;
        d->nstate  = 2;
    }
    else
    {
        d->nparams = 4*(1+ninter);
        d->nstate  = 2+ninter;
    }
    d->bSum      = bSum;
    d->bDiscrete = bDiscrete;
    snew(d->beta, d->nreplica);
    snew(d->state, d->nreplica);
    snew(d->bMask, d->nreplica);
    snew(d->d2_replica, d->nreplica);
    snew(d->sumft, d->nframe);
    snew(d->sumit, d->nframe);
    snew(d->sumfct, d->nframe);
    snew(d->sumict, d->nframe);
    snew(d->params, d->nparams);
    snew(d->fcalt, d->nreplica);
    snew(d->icalt, d->nreplica);

    /* convert_times(d->nframe,d->time); */

    if (t0 < 0)
    {
        d->j0 = 0;
    }
    else
    {
        for (d->j0 = 0; (d->j0 < d->nframe) && (d->time[d->j0] < t0); d->j0++)
        {
            ;
        }
    }
    if (t1 < 0)
    {
        d->j1 = d->nframe;
    }
    else
    {
        for (d->j1 = 0; (d->j1 < d->nframe) && (d->time[d->j1] < t1); d->j1++)
        {
            ;
        }
    }
    if ((d->j1-d->j0) < d->nparams+2)
    {
        gmx_fatal(FARGS, "Start (%f) or end time (%f) for fitting inconsistent. Reduce t0, increase t1 or supply more data", t0, t1);
    }
    fprintf(fp, "Will optimize from %g to %g\n",
            d->time[d->j0], d->time[d->j1-1]);
    d->nmask = d->nreplica;
    for (i = 0; (i < d->nreplica); i++)
    {
        snew(d->beta[i], d->nframe);
        snew(d->state[i], d->nframe);
        snew(d->fcalt[i], d->nframe);
        snew(d->icalt[i], d->nframe);
        d->bMask[i] = TRUE;
        for (j = 0; (j < d->nframe); j++)
        {
            d->beta[i][j] = 1.0/(BOLTZ*d->temp[i][j]);
            dd            = d->data[i][j];
            if (bDiscrete)
            {
                if (dd <= cutoff)
                {
                    d->state[i][j] = 0;
                }
                else if ((ucut > cutoff) && (dd <= ucut))
                {
                    d->state[i][j] = 1;
                }
                else
                {
                    d->state[i][j] = d->nstate-1;
                }
            }
            else
            {
                d->state[i][j] = dd*nmult;
            }
        }
    }
    sum_ft(d);

    /* Assume forward rate constant is half the total time in this
     * simulation and backward is ten times as long */
    if (bDiscrete)
    {
        tau_f            = d->time[d->nframe-1];
        tau_u            = 4*tau_f;
        d->params[epEuf] = Euf;
        d->params[epAuf] = exp(d->params[epEuf]/(BOLTZ*tref))/tau_f;
        if (bBack)
        {
            d->params[epEfu] = Efu;
            d->params[epAfu] = exp(d->params[epEfu]/(BOLTZ*tref))/tau_u;
            if (ninter > 0)
            {
                d->params[eqEui] = Ei;
                d->params[eqAui] = exp(d->params[eqEui]/(BOLTZ*tref))/tau_u;
                d->params[eqEiu] = Ei;
                d->params[eqAiu] = exp(d->params[eqEiu]/(BOLTZ*tref))/tau_u;
            }
        }
        else
        {
            d->params[epAfu]  = 0;
            d->params[epEfu]  = 0;
        }
    }
    else
    {
        d->params[epEuf] = Euf;
        if (d->data2)
        {
            d->params[epAuf] = 0.1;
        }
        else
        {
            d->params[epAuf] = 20.0;
        }
    }
}

static real tau(real A, real E, real T)
{
    return exp(E/(BOLTZ*T))/A;
}

static real folded_fraction(t_remd_data *d, real tref)
{
    real tauf, taub;

    tauf = tau(d->params[epAuf], d->params[epEuf], tref);
    taub = tau(d->params[epAfu], d->params[epEfu], tref);

    return (taub/(tauf+taub));
}

static void print_tau(FILE *gp, t_remd_data *d, real tref)
{
    real tauf, taub, ddd, fff, DG, DH, TDS, Tm, Tb, Te, Fb, Fe, Fm;
    int  i, np = d->nparams;

    ddd = calc_d2(d);
    fprintf(gp, "Final value for Chi2 = %12.5e (%d replicas)\n", ddd, d->nmask);
    tauf = tau(d->params[epAuf], d->params[epEuf], tref);
    fprintf(gp, "%s = %12.5e %s = %12.5e (kJ/mole)\n",
            epnm(np, epAuf), d->params[epAuf],
            epnm(np, epEuf), d->params[epEuf]);
    if (bBack(d))
    {
        taub = tau(d->params[epAfu], d->params[epEfu], tref);
        fprintf(gp, "%s = %12.5e %s = %12.5e (kJ/mole)\n",
                epnm(np, epAfu), d->params[epAfu],
                epnm(np, epEfu), d->params[epEfu]);
        fprintf(gp, "Equilibrium properties at T = %g\n", tref);
        fprintf(gp, "tau_f = %8.3f ns, tau_b = %8.3f ns\n", tauf/1000, taub/1000);
        fff = taub/(tauf+taub);
        DG  = BOLTZ*tref*log(fff/(1-fff));
        DH  = d->params[epEfu]-d->params[epEuf];
        TDS = DH-DG;
        fprintf(gp, "Folded fraction     F = %8.3f\n", fff);
        fprintf(gp, "Unfolding energies: DG = %8.3f  DH = %8.3f TDS = %8.3f\n",
                DG, DH, TDS);
        Tb = 260;
        Te = 420;
        Tm = 0;
        Fm = 0;
        Fb = folded_fraction(d, Tb);
        Fe = folded_fraction(d, Te);
        while ((Te-Tb > 0.001) && (Fm != 0.5))
        {
            Tm = 0.5*(Tb+Te);
            Fm = folded_fraction(d, Tm);
            if (Fm > 0.5)
            {
                Fb = Fm;
                Tb = Tm;
            }
            else if (Fm < 0.5)
            {
                Te = Tm;
                Fe = Fm;
            }
        }
        if ((Fb-0.5)*(Fe-0.5) <= 0)
        {
            fprintf(gp, "Melting temperature Tm = %8.3f K\n", Tm);
        }
        else
        {
            fprintf(gp, "No melting temperature detected between 260 and 420K\n");
        }
        if (np > 4)
        {
            char *ptr;
            fprintf(gp, "Data for intermediates at T = %g\n", tref);
            fprintf(gp, "%8s  %10s  %10s  %10s\n", "Name", "A", "E", "tau");
            for (i = 0; (i < np/2); i++)
            {
                tauf = tau(d->params[2*i], d->params[2*i+1], tref);
                ptr  = epnm(d->nparams, 2*i);
                fprintf(gp, "%8s  %10.3e  %10.3e  %10.3e\n", ptr+1,
                        d->params[2*i], d->params[2*i+1], tauf/1000);
            }
        }
    }
    else
    {
        fprintf(gp, "Equilibrium properties at T = %g\n", tref);
        fprintf(gp, "tau_f = %8.3f\n", tauf);
    }
}

static void dump_remd_parameters(FILE *gp, t_remd_data *d, const char *fn,
                                 const char *fn2, const char *rfn,
                                 const char *efn, const char *mfn, int skip, real tref,
                                 output_env_t oenv)
{
    FILE       *fp, *hp;
    int         i, j, np = d->nparams;
    real        rhs, tauf, taub, fff, DG;
    real       *params;
    const char *leg[]  = { "Measured", "Fit", "Difference" };
    const char *mleg[] = { "Folded fraction", "DG (kJ/mole)"};
    char      **rleg;
    real        fac[] = { 0.97, 0.98, 0.99, 1.0, 1.01, 1.02, 1.03 };
#define NFAC asize(fac)
    real        d2[NFAC];
    double      norm;

    integrate_dfdt(d);
    print_tau(gp, d, tref);
    norm = (d->bDiscrete ? 1.0/d->nmask : 1.0);

    if (fn)
    {
        fp = xvgropen(fn, "Optimized fit to data", "Time (ps)", "Fraction Folded", oenv);
        xvgr_legend(fp, asize(leg), leg, oenv);
        for (i = 0; (i < d->nframe); i++)
        {
            if ((skip <= 0) || ((i % skip) == 0))
            {
                fprintf(fp, "%12.5e  %12.5e  %12.5e  %12.5e\n", d->time[i],
                        d->sumft[i]*norm, d->sumfct[i]*norm,
                        (d->sumft[i]-d->sumfct[i])*norm);
            }
        }
        ffclose(fp);
    }
    if (!d->bSum && rfn)
    {
        snew(rleg, d->nreplica*2);
        for (i = 0; (i < d->nreplica); i++)
        {
            snew(rleg[2*i], 32);
            snew(rleg[2*i+1], 32);
            sprintf(rleg[2*i], "\\f{4}F(t) %d", i);
            sprintf(rleg[2*i+1], "\\f{12}F \\f{4}(t) %d", i);
        }
        fp = xvgropen(rfn, "Optimized fit to data", "Time (ps)", "Fraction Folded", oenv);
        xvgr_legend(fp, d->nreplica*2, (const char**)rleg, oenv);
        for (j = 0; (j < d->nframe); j++)
        {
            if ((skip <= 0) || ((j % skip) == 0))
            {
                fprintf(fp, "%12.5e", d->time[j]);
                for (i = 0; (i < d->nreplica); i++)
                {
                    fprintf(fp, "  %5f  %9.2e", is_folded(d, i, j), d->fcalt[i][j]);
                }
                fprintf(fp, "\n");
            }
        }
        ffclose(fp);
    }

    if (fn2 && (d->nstate > 2))
    {
        fp = xvgropen(fn2, "Optimized fit to data", "Time (ps)",
                      "Fraction Intermediate", oenv);
        xvgr_legend(fp, asize(leg), leg, oenv);
        for (i = 0; (i < d->nframe); i++)
        {
            if ((skip <= 0) || ((i % skip) == 0))
            {
                fprintf(fp, "%12.5e  %12.5e  %12.5e  %12.5e\n", d->time[i],
                        d->sumit[i]*norm, d->sumict[i]*norm,
                        (d->sumit[i]-d->sumict[i])*norm);
            }
        }
        ffclose(fp);
    }
    if (mfn)
    {
        if (bBack(d))
        {
            fp = xvgropen(mfn, "Melting curve", "T (K)", "", oenv);
            xvgr_legend(fp, asize(mleg), mleg, oenv);
            for (i = 260; (i <= 420); i++)
            {
                tauf = tau(d->params[epAuf], d->params[epEuf], 1.0*i);
                taub = tau(d->params[epAfu], d->params[epEfu], 1.0*i);
                fff  = taub/(tauf+taub);
                DG   = BOLTZ*i*log(fff/(1-fff));
                fprintf(fp, "%5d  %8.3f  %8.3f\n", i, fff, DG);
            }
            ffclose(fp);
        }
    }

    if (efn)
    {
        snew(params, d->nparams);
        for (i = 0; (i < d->nparams); i++)
        {
            params[i] = d->params[i];
        }

        hp = xvgropen(efn, "Chi2 as a function of relative parameter",
                      "Fraction", "Chi2", oenv);
        for (j = 0; (j < d->nparams); j++)
        {
            /* Reset all parameters to optimized values */
            fprintf(hp, "@type xy\n");
            for (i = 0; (i < d->nparams); i++)
            {
                d->params[i] = params[i];
            }
            /* Now modify one of them */
            for (i = 0; (i < NFAC); i++)
            {
                d->params[j] = fac[i]*params[j];
                d2[i]        = calc_d2(d);
                fprintf(gp, "%s = %12g  d2 = %12g\n", epnm(np, j), d->params[j], d2[i]);
                fprintf(hp, "%12g  %12g\n", fac[i], d2[i]);
            }
            fprintf(hp, "&\n");
        }
        ffclose(hp);
        for (i = 0; (i < d->nparams); i++)
        {
            d->params[i] = params[i];
        }
        sfree(params);
    }
    if (!d->bSum)
    {
        for (i = 0; (i < d->nreplica); i++)
        {
            fprintf(gp, "Chi2[%3d] = %8.2e\n", i, d->d2_replica[i]);
        }
    }
}
#endif /*HAVE_LIBGSL*/

int gmx_kinetics(int argc, char *argv[])
{
    const char     *desc[] = {
        "[TT]g_kinetics[tt] reads two [TT].xvg[tt] files, each one containing data for N replicas.",
        "The first file contains the temperature of each replica at each timestep,",
        "and the second contains real values that can be interpreted as",
        "an indicator for folding. If the value in the file is larger than",
        "the cutoff it is taken to be unfolded and the other way around.[PAR]",
        "From these data an estimate of the forward and backward rate constants",
        "for folding is made at a reference temperature. In addition,",
        "a theoretical melting curve and free energy as a function of temperature",
        "are printed in an [TT].xvg[tt] file.[PAR]",
        "The user can give a max value to be regarded as intermediate",
        "([TT]-ucut[tt]), which, when given will trigger the use of an intermediate state",
        "in the algorithm to be defined as those structures that have",
        "cutoff < DATA < ucut. Structures with DATA values larger than ucut will",
        "not be regarded as potential folders. In this case 8 parameters are optimized.[PAR]",
        "The average fraction foled is printed in an [TT].xvg[tt] file together with the fit to it.",
        "If an intermediate is used a further file will show the build of the intermediate and the fit to that process.[PAR]",
        "The program can also be used with continuous variables (by setting",
        "[TT]-nodiscrete[tt]). In this case kinetics of other processes can be",
        "studied. This is very much a work in progress and hence the manual",
        "(this information) is lagging behind somewhat.[PAR]",
        "In order to compile this program you need access to the GNU",
        "scientific library."
    };
    static int      nreplica  = 1;
    static real     tref      = 298.15;
    static real     cutoff    = 0.2;
    static real     ucut      = 0.0;
    static real     Euf       = 10;
    static real     Efu       = 30;
    static real     Ei        = 10;
    static gmx_bool bHaveT    = TRUE;
    static real     t0        = -1;
    static real     t1        = -1;
    static real     tb        = 0;
    static real     te        = 0;
    static real     tol       = 1e-3;
    static int      maxiter   = 100;
    static int      skip      = 0;
    static int      nmult     = 1;
    static gmx_bool bBack     = TRUE;
    static gmx_bool bSplit    = TRUE;
    static gmx_bool bSum      = TRUE;
    static gmx_bool bDiscrete = TRUE;
    t_pargs         pa[]      = {
        { "-time",    FALSE, etBOOL, {&bHaveT},
          "Expect a time in the input" },
        { "-b",       FALSE, etREAL, {&tb},
          "First time to read from set" },
        { "-e",       FALSE, etREAL, {&te},
          "Last time to read from set" },
        { "-bfit",    FALSE, etREAL, {&t0},
          "Time to start the fit from" },
        { "-efit",    FALSE, etREAL, {&t1},
          "Time to end the fit" },
        { "-T",       FALSE, etREAL, {&tref},
          "Reference temperature for computing rate constants" },
        { "-n",       FALSE, etINT, {&nreplica},
          "Read data for this number of replicas. Only necessary when files are written in xmgrace format using @type and & as delimiters." },
        { "-cut",     FALSE, etREAL, {&cutoff},
          "Cut-off (max) value for regarding a structure as folded" },
        { "-ucut",    FALSE, etREAL, {&ucut},
          "Cut-off (max) value for regarding a structure as intermediate (if not folded)" },
        { "-euf",     FALSE, etREAL, {&Euf},
          "Initial guess for energy of activation for folding (kJ/mol)" },
        { "-efu",     FALSE, etREAL, {&Efu},
          "Initial guess for energy of activation for unfolding (kJ/mol)" },
        { "-ei",      FALSE, etREAL, {&Ei},
          "Initial guess for energy of activation for intermediates (kJ/mol)" },
        { "-maxiter", FALSE, etINT, {&maxiter},
          "Max number of iterations" },
        { "-back",    FALSE, etBOOL, {&bBack},
          "Take the back reaction into account" },
        { "-tol",     FALSE, etREAL, {&tol},
          "Absolute tolerance for convergence of the Nelder and Mead simplex algorithm" },
        { "-skip",    FALSE, etINT, {&skip},
          "Skip points in the output [TT].xvg[tt] file" },
        { "-split",   FALSE, etBOOL, {&bSplit},
          "Estimate error by splitting the number of replicas in two and refitting" },
        { "-sum",     FALSE, etBOOL, {&bSum},
          "Average folding before computing [GRK]chi[grk]^2" },
        { "-discrete", FALSE, etBOOL, {&bDiscrete},
          "Use a discrete folding criterion (F <-> U) or a continuous one" },
        { "-mult",    FALSE, etINT, {&nmult},
          "Factor to multiply the data with before discretization" }
    };
#define NPA asize(pa)

    FILE        *fp;
    real         dt_t, dt_d, dt_d2;
    int          nset_t, nset_d, nset_d2, n_t, n_d, n_d2, i;
    const char  *tfile, *dfile, *dfile2;
    t_remd_data  remd;
    output_env_t oenv;

    t_filenm     fnm[] = {
        { efXVG, "-f",    "temp",    ffREAD   },
        { efXVG, "-d",    "data",    ffREAD   },
        { efXVG, "-d2",   "data2",   ffOPTRD  },
        { efXVG, "-o",    "ft_all",  ffWRITE  },
        { efXVG, "-o2",   "it_all",  ffOPTWR  },
        { efXVG, "-o3",   "ft_repl", ffOPTWR  },
        { efXVG, "-ee",   "err_est", ffOPTWR  },
        { efLOG, "-g",    "remd",    ffWRITE  },
        { efXVG, "-m",    "melt",    ffWRITE  }
    };
#define NFILE asize(fnm)

    CopyRight(stderr, argv[0]);
    parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_BE_NICE | PCA_TIME_UNIT,
                      NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv);

#ifdef HAVE_LIBGSL
    please_cite(stdout, "Spoel2006d");
    if (cutoff < 0)
    {
        gmx_fatal(FARGS, "cutoff should be >= 0 (rather than %f)", cutoff);
    }

    tfile   = opt2fn("-f", NFILE, fnm);
    dfile   = opt2fn("-d", NFILE, fnm);
    dfile2  = opt2fn_null("-d2", NFILE, fnm);

    fp = ffopen(opt2fn("-g", NFILE, fnm), "w");

    remd.temp = read_xvg_time(tfile, bHaveT,
                              opt2parg_bSet("-b", NPA, pa), tb,
                              opt2parg_bSet("-e", NPA, pa), te,
                              nreplica, &nset_t, &n_t, &dt_t, &remd.time);
    printf("Read %d sets of %d points in %s, dt = %g\n\n", nset_t, n_t, tfile, dt_t);
    sfree(remd.time);

    remd.data = read_xvg_time(dfile, bHaveT,
                              opt2parg_bSet("-b", NPA, pa), tb,
                              opt2parg_bSet("-e", NPA, pa), te,
                              nreplica, &nset_d, &n_d, &dt_d, &remd.time);
    printf("Read %d sets of %d points in %s, dt = %g\n\n", nset_d, n_d, dfile, dt_d);

    if ((nset_t != nset_d) || (n_t != n_d) || (dt_t != dt_d))
    {
        gmx_fatal(FARGS, "Files %s and %s are inconsistent. Check log file",
                  tfile, dfile);
    }

    if (dfile2)
    {
        remd.data2 = read_xvg_time(dfile2, bHaveT,
                                   opt2parg_bSet("-b", NPA, pa), tb,
                                   opt2parg_bSet("-e", NPA, pa), te,
                                   nreplica, &nset_d2, &n_d2, &dt_d2, &remd.time);
        printf("Read %d sets of %d points in %s, dt = %g\n\n",
               nset_d2, n_d2, dfile2, dt_d2);
        if ((nset_d2 != nset_d) || (n_d != n_d2) || (dt_d != dt_d2))
        {
            gmx_fatal(FARGS, "Files %s and %s are inconsistent. Check log file",
                      dfile, dfile2);
        }
    }
    else
    {
        remd.data2 = NULL;
    }

    remd.nreplica  = nset_d;
    remd.nframe    = n_d;
    remd.dt        = 1;
    preprocess_remd(fp, &remd, cutoff, tref, ucut, bBack, Euf, Efu, Ei, t0, t1,
                    bSum, bDiscrete, nmult);

    optimize_remd_parameters(fp, &remd, maxiter, tol);

    dump_remd_parameters(fp, &remd, opt2fn("-o", NFILE, fnm),
                         opt2fn_null("-o2", NFILE, fnm),
                         opt2fn_null("-o3", NFILE, fnm),
                         opt2fn_null("-ee", NFILE, fnm),
                         opt2fn("-m", NFILE, fnm), skip, tref, oenv);

    if (bSplit)
    {
        printf("Splitting set of replicas in two halves\n");
        for (i = 0; (i < remd.nreplica); i++)
        {
            remd.bMask[i] = FALSE;
        }
        remd.nmask = 0;
        for (i = 0; (i < remd.nreplica); i += 2)
        {
            remd.bMask[i] = TRUE;
            remd.nmask++;
        }
        sum_ft(&remd);
        optimize_remd_parameters(fp, &remd, maxiter, tol);
        dump_remd_parameters(fp, &remd, "test1.xvg", NULL, NULL, NULL, NULL, skip, tref, oenv);

        for (i = 0; (i < remd.nreplica); i++)
        {
            remd.bMask[i] = !remd.bMask[i];
        }
        remd.nmask = remd.nreplica - remd.nmask;

        sum_ft(&remd);
        optimize_remd_parameters(fp, &remd, maxiter, tol);
        dump_remd_parameters(fp, &remd, "test2.xvg", NULL, NULL, NULL, NULL, skip, tref, oenv);

        for (i = 0; (i < remd.nreplica); i++)
        {
            remd.bMask[i] = FALSE;
        }
        remd.nmask = 0;
        for (i = 0; (i < remd.nreplica/2); i++)
        {
            remd.bMask[i] = TRUE;
            remd.nmask++;
        }
        sum_ft(&remd);
        optimize_remd_parameters(fp, &remd, maxiter, tol);
        dump_remd_parameters(fp, &remd, "test1.xvg", NULL, NULL, NULL, NULL, skip, tref, oenv);

        for (i = 0; (i < remd.nreplica); i++)
        {
            remd.bMask[i] = FALSE;
        }
        remd.nmask = 0;
        for (i = remd.nreplica/2; (i < remd.nreplica); i++)
        {
            remd.bMask[i] = TRUE;
            remd.nmask++;
        }
        sum_ft(&remd);
        optimize_remd_parameters(fp, &remd, maxiter, tol);
        dump_remd_parameters(fp, &remd, "test1.xvg", NULL, NULL, NULL, NULL, skip, tref, oenv);
    }
    ffclose(fp);

    view_all(oenv, NFILE, fnm);

    thanx(stderr);
#else
    fprintf(stderr, "This program should be compiled with the GNU scientific library. Please install the library and reinstall GROMACS.\n");
#endif /*HAVE_LIBGSL*/

    return 0;
}


text:
/*
** conversions: p->z, z->p, t->z
**
** G.Lohmann, April 2007
*/
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <stdlib.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_errno.h>

#define ABS(x) ((x) > 0 ? (x) : -(x))

/*
** convert t to p values
*/
double t2p(double t,double df)
{
  double a,b,x;
  extern double gsl_sf_beta_inc(double,double,double);

  x = df/(df+(t*t));
  if (x < 0 || x > 1) return 1;
  a = 0.5*df;
  b = 0.5;

  return 0.5 * gsl_sf_beta_inc(a,b,x);
}

/*
** convert t to z values
*/
double t2z(double t,double df)
{
  double p=0,z=0;
  double a,b,x;
  extern double gsl_sf_beta_inc(double,double,double);

  /* t to p */
  x = df/(df+(t*t));
  if (x <= 0 || x > 1) return 0;

  a = 0.5*df;
  b = 0.5;
  p = 0.5 * gsl_sf_beta_inc(a,b,x);

  /* p to z */
  z = gsl_cdf_ugaussian_Qinv(p);
  return z;
}

/*
** approximation: convert t to p values
*/
float
t2z_approx(float t,float df)
{
  float z=0,u;

  u = df*log(1.0+t*t/df)*(1.0-0.5/df);
  if (u <= 0) return 0;
  z = sqrt(u);
  return z;
}

/*
** convert p to t values
*/
double
p2t(double px,double df)
{
  double p,t0,t1,t,step=0.00001;

  t  = 0;
  t0 = 0;
  t1 = 20;
  while (ABS(t0-t1) > step) {
    t = (t0 + t1)*0.5;
    p = t2p(t,df);
    if (p < px)  t1 = t;
    else t0 = t;
  }
  p = t2p(t,df);

  return t;
}

/*
** convert z to p values
*/
double z2p(double z)
{
  if (z < 0)
    return gsl_cdf_ugaussian_Q(-z);
  else
    return gsl_cdf_ugaussian_Q(z);
}

/*
** convert p to z values
*/
double p2z(double p)
{
  return gsl_cdf_ugaussian_Qinv(p);
}

ext: c

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
/** \file real-compressor.c */
#ifndef __COMPLEARN_REALCOMPRESSOR_H
#define __COMPLEARN_REALCOMPRESSOR_H

#define COMPLEARN_REAL_COMPRESSOR_TYPE        (real_compressor_get_type ())
#define COMPLEARN_TYPE_REAL_COMPRESSOR        (real_compressor_get_type ())
#define COMPLEARN_REAL_COMPRESSOR(obj)        (G_TYPE_CHECK_INSTANCE_CAST ((obj), COMPLEARN_REAL_COMPRESSOR_TYPE, CompLearnRealCompressor))
#define IS_COMPLEARN_REAL_COMPRESSOR(obj)        (G_TYPE_CHECK_INSTANCE_TYPE ((obj), COMPLEARN_REAL_COMPRESSOR_TYPE))
#define COMPLEARN_REAL_COMPRESSOR_GET_INTERFACE(inst) (G_TYPE_INSTANCE_GET_INTERFACE((inst), COMPLEARN_REAL_COMPRESSOR_TYPE, CompLearnRealCompressorIface))

#include <glib.h>
#include <glib-object.h>
#include <glib/garray.h>
#include <gsl/gsl_matrix.h>

#define COMPLEARN_ERROR 1
#define COMPLEARN_ERROR_NO_COMPRESSOR_SET 100

typedef struct _CompLearnRealCompressor CompLearnRealCompressor;
typedef struct _CompLearnRealCompressorIface CompLearnRealCompressorIface;

struct _CompLearnRealCompressorIface {
  GTypeInterface parent;

  GString *(*compress)(CompLearnRealCompressor *self, const GString *input);
  GString *(*decompress)(CompLearnRealCompressor *self, const GString *input);
  GString *(*blurb)(CompLearnRealCompressor *self);
  GString *(*canonical_extension)(CompLearnRealCompressor *self);
  GString *(*name)(CompLearnRealCompressor *self);
  GString *(*compressor_version)(CompLearnRealCompressor *self);
  GString *(*binding_version)(CompLearnRealCompressor *self);
  gboolean (*is_threadsafe)(CompLearnRealCompressor *self);
  gboolean (*is_compressible)(CompLearnRealCompressor *self, const GString *input);
  gboolean (*is_decompressible)(CompLearnRealCompressor *self, const GString *input);
  gboolean (*is_just_size)(CompLearnRealCompressor *self);
  gboolean (*is_hash_function)(CompLearnRealCompressor *self);
  GString *(*hash)(CompLearnRealCompressor *self, const GString *input);
  gboolean (*is_operational)(CompLearnRealCompressor *self);
  gboolean (*is_private_property)(CompLearnRealCompressor *self, const char *propname);
  gdouble (*compressed_size)(CompLearnRealCompressor *self, const GString *input);
  guint64 (*window_size)(CompLearnRealCompressor *self);
  CompLearnRealCompressor *(*clone)(CompLearnRealCompressor *self);
};

GType real_compressor_get_type(void);

/// Compress the string by the compressor
/// \param[in] self is a pointer to compressor
/// \param[in] input is a pointer to a GString which must be compressed
/// \return a pointer to compressed string
GString *real_compressor_compress(CompLearnRealCompressor *self,const GString *input);
GString *real_compressor_hash(CompLearnRealCompressor *self,const GString *input);
/// Decompress the string by the compressor
/// \param[in] self is pointer to compressor
/// \param[in] input is a pointer to a GString which must be decompressed
/// \return a pointer to decompressed string
GString *real_compressor_decompress(CompLearnRealCompressor *self,const GString *input);
/// Calculate size of the compressed string
/// \param[in] self is a pointer to a compressor
/// \param[in] input is a pointer to string which must be compress
/// \return size of compressed string in bits
gdouble real_compressor_compressed_size(CompLearnRealCompressor *self,const GString *input);
/// Return a short text description of this compressor
/// \param[in] self is a pointer to a compressor
GString *real_compressor_blurb(CompLearnRealCompressor *self);
GString *real_compressor_name(CompLearnRealCompressor *self);
/// Display the compressor version
/// param[in] is a pointer to compressor
/// \return a GString holding a 3-part version
GString *real_compressor_compressor_version(CompLearnRealCompressor *self);
/// Display the compressor binding version
/// \param[in] is a pointer to compressor
/// \return a GString holding a 3-part version
GString *real_compressor_binding_version(CompLearnRealCompressor *self);
/// This function checks whether the string is compressible by the compressor
/// \param[in] self is a pointer to compressor
/// \param[in] input is a pointer to GString which is under test
/// \return true for every string
gboolean real_compressor_is_compressible(CompLearnRealCompressor *self, const GString *input);
/// Check whether the string is decompressible by the compressor
/// \param[in] self is a pointer to compressor
/// \param[in] input is a pointer to string which is under test
gboolean real_compressor_is_decompressible(CompLearnRealCompressor *self, const GString *input);

gboolean real_compressor_is_private_property(CompLearnRealCompressor *self, const char *input);
guint64 real_compressor_window_size(CompLearnRealCompressor *self);
gboolean real_compressor_is_threadsafe(CompLearnRealCompressor *self);
gboolean real_compressor_is_just_size(CompLearnRealCompressor *self);
gboolean real_compressor_is_hash_function(CompLearnRealCompressor *self);
GString *real_compressor_canonical_extension(CompLearnRealCompressor *rc);
gboolean real_compressor_is_operational(CompLearnRealCompressor *self);
CompLearnRealCompressor *real_compressor_clone(CompLearnRealCompressor *self);

#define SET_DEFAULT_PROPS(groupname, clt, mobj) \
  do { \
  GParamSpec **gps, **cur; \
  g_assert(mobj != NULL); \
    if (complearn_environment_get_nameable(groupname) == NULL) { \
      complearn_environment_register_nameable(groupname, G_OBJECT(mobj)); \
      break; \
    } \
  gps =g_object_class_list_properties(G_OBJECT_CLASS(clt(mobj)), NULL); \
  for (cur = gps; *cur; cur += 1) { \
    GValue v = {0,}; \
    g_value_init(&v, (*cur)->value_type); \
    g_param_value_set_default(*cur, &v); \
    g_object_set_property(G_OBJECT(mobj), (*cur)->name, &v); \
    complearn_environment_register_property(G_OBJECT(mobj), (*cur), &v); \
  } } while(0)

#define SET_DEFAULT_COMPRESSOR_PROPS(groupname, clt, mobj) \
  do { \
      SET_DEFAULT_PROPS(groupname, clt, mobj); \
    if (complearn_environment_get_nameable(groupname) == NULL) { \
      complearn_environment_register_compressor(mobj); \
    } \
  } while (0);

#endif

#define G_LOG_LEVEL_NOTICE G_LOG_LEVEL_USER_SHIFT
#define g_notice(...) g_log(G_LOG_DOMAIN, G_LOG_LEVEL_NOTICE, __VA_ARGS__)

void real_compressor_interface_init (gpointer g_iface, gpointer iface_data);

max_stars_repo_name: rudi-cilibrasi/classic-complearn

max_issues_repo_name: rudi-cilibrasi/classic-complearn

max_forks_repo_name: rudi-cilibrasi/classic-complearn


text:
/*
** Read design file
**
** G.Lohmann, MPI-KYB 2018
*/
#include <viaio/Vlib.h>
#include <viaio/VImage.h>
#include <viaio/mu.h>
#include <viaio/option.h>

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>

typedef struct TrialStruct {
  int   id;
  float onset;
  float duration;
  float height;
} Trial;

extern void printmat(gsl_matrix *R,char *str);
extern void printvec(gsl_vector *x,char *str);

int test_ascii(int val)
{
  if (val >= 'a' && val <= 'z') return 1;
  if (val >= 'A' && val <= 'Z') return 1;
  if (val >= '0' && val <= '9') return 1;
  if (val ==  ' ') return 1;
  if (val ==  '-') return 1;
  if (val ==  '+') return 1;
  if (val ==  '.') return 1;
  if (val == '\0') return 1;
  if (val == '\n') return 1;
  if (val == '\r') return 1;
  if (val == '\t') return 1;
  if (val == '\v') return 1;
  return 0;
}

int line_empty(char *buf,int len)
{
  int i;
  for (i=0; i<len; i++) {
    if (buf[i] != ' ' && buf[i] != '\0' && buf[i] != '\n') return 0;
  }
  return 1;
}

int CheckBuffer(char *buf,int len)
{
  int j;
  if(strlen(buf) < 1) return 0;
  if (buf[0] == '%' || buf[0] == '#' || buf[0] == '/') return 0;  /* comment */
  for (j=0; j<len; j++) {
    if (buf[j] == '\v') buf[j] = ' '; /* remove tabs */
    if (buf[j] == '\t') buf[j] = ' ';
  }
  if (line_empty(buf,len) > 0) return 0;
  return 1;
}

int VistaFormat(char *buf,int len)
{
  if (strncmp(buf,"V-data",6) == 0) return 1;
  return 0;
}

/* parse design file */
Trial *ReadDesign(VString designfile,int *numtrials,int *nevents)
{
  FILE *fp=NULL;
  int  i,j,id,len=4096;
  float onset=0,duration=0,height=0;
  char *buf = (char *)VCalloc(len,sizeof(char));

  fp = fopen(designfile,"r");
  if (!fp) VError(" error opening design file %s",designfile);
  int ntrials=0;
  while (!feof(fp)) {
    for (j=0; j<len; j++) buf[j] = '\0';
    if (fgets(buf,len,fp) == NULL) break;
    if (VistaFormat(buf,len) > 1) VError(" Design file must be a text file");
    if (CheckBuffer(buf,len) < 1) continue;
    if (! test_ascii((int)buf[0])) VError(" File %s: line %d begins with an illegal character (%c)",designfile,ntrials,buf[0]);
    ntrials++;
  }
  *numtrials = ntrials;

  Trial *trial = (Trial *) VCalloc(ntrials+1,sizeof(Trial));

  i = (*nevents) = 0;
  rewind(fp);
  while (!feof(fp)) {
    for (j=0; j<len; j++) buf[j] = '\0';
    if (fgets(buf,len,fp) == NULL) break;
    if (CheckBuffer(buf,len) < 1) continue;

    int nch = sscanf(buf,"%d %f %f %f",&id,&onset,&duration,&height);
    if (nch > 0 && nch != 4) VError(" line %d: illegal input format",i+1);
    if (duration < 0.5 && duration >= -0.0001) duration = 0.5;
    trial[i].id       = id;
    trial[i].onset    = onset;
    trial[i].duration = duration;
    trial[i].height   = height;
    i++;
    if (id > (*nevents)) (*nevents) = id;
  }
  fclose(fp);
  return trial;
}

Trial *CopyTrials(Trial *trial,int numtrials)
{
  int i;
  Trial *newtrial = (Trial *) VCalloc(numtrials,sizeof(Trial));
  for (i=0; i<numtrials; i++) {
    newtrial[i].id       = trial[i].id;
    newtrial[i].onset    = trial[i].onset;
    newtrial[i].duration = trial[i].duration;
    newtrial[i].height   = trial[i].height;
  }
  return newtrial;
}

/* concatenate trials from all runs */
Trial *ConcatenateTrials(Trial **trial,int *numtrials,float *run_duration,int dlists,int sumtrials)
{
  int i,j;
  Trial *alltrials = (Trial *) VCalloc(sumtrials,sizeof(Trial));

  int ii=0;
  float add=0;
  for (i=0; i<dlists; i++) {
    for (j=0; j<numtrials[i]; j++) {
      alltrials[ii].id       = trial[i][j].id;
      alltrials[ii].onset    = trial[i][j].onset + add;
      alltrials[ii].duration = trial[i][j].duration;
      alltrials[ii].height   = trial[i][j].height;
      ii++;
    }
    add += run_duration[i];
  }
  return alltrials;
}

/* check if trial labels are complete and positive */
void CheckTrialLabels(Trial *trial,int numtrials)
{
  int i,j;
  int maxlabel = -1;
  for (i=0; i<numtrials; i++) {
    j = trial[i].id;
    if (j > maxlabel) maxlabel = j;
    if (j < 1) VError(" trial[%d] has illegal label %d",i,j);
  }

  int *tmp = (int *) VCalloc(maxlabel+1,sizeof(int));
  for (i=0; i<numtrials; i++) {
    j = trial[i].id;
    tmp[j]++;
  }
  for (i=1; i<=maxlabel; i++) {
    if (tmp[i] < 1) VError(" Label %d missing in design file",i);
    if (tmp[i] < 4) VWarning(" Label %d has too few trials (%d) for random permutations",i,tmp[i]);
  }
  VFree(tmp);
}


text:
/*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*
**                                                                                **
** This file forms part of the Underworld geophysics modelling application.         **
**                                                                                **
** For full license and copyright information, please refer to the LICENSE.md file  **
** located at the project root, or contact the authors.                             **
**                                                                                **
**~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*/
#include <petsc.h>
#include <petscvec.h>
#include <petscmat.h>
#include <petscksp.h>

#include "common-driver-utils.h"

/*

writes out

file1
  1) its, residual, time
  2) a summary of results
file2
  1) solver configuration

*/
PetscErrorCode BSSCR_solver_output( KSP ksp, PetscInt monitor_index, const char res_file[], const char solver_conf[] )
{
	PetscViewer v;
	MPI_Comm comm;
	PetscInt i,nargs;
	char **args;

	PetscObjectGetComm( (PetscObject)ksp, &comm );
	PetscViewerASCIIOpen( comm, res_file, &v );

	BSSCR_GeneratePetscHeader_for_viewer( v );
        PetscViewerASCIIPrintf( v, "\n");
	BSSCR_KSPLogSolve( v, monitor_index, ksp );
	BSSCR_BSSCR_KSPLogSolveSummary( v, monitor_index, ksp );

	Stg_PetscViewerDestroy(&v );

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

	PetscObjectGetComm( (PetscObject)ksp, &comm );
	PetscViewerASCIIOpen( comm, solver_conf, &v );

	BSSCR_GeneratePetscHeader_for_viewer( v );
	PetscViewerASCIIPrintf( v, "\n\nSolver configuration: \n");

	PetscViewerASCIIPrintf( v, "\n1) Command line arguments supplied: \n");
	PetscViewerASCIIPushTab(v);
        PetscGetArgs( &nargs, &args );
        for( i=0; i<nargs; i++ ) {
                PetscViewerASCIIPrintf( v, "[%d]: %s \n", i, args[i] );
        }
	PetscViewerASCIIPopTab(v);

	PetscViewerASCIIPrintf( v, "\n2) Petsc KSPView(): \n");
	KSPView( ksp, v );

	Stg_PetscViewerDestroy(&v );

	PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "BSSCR_KSPLogConvergenceRate"
PetscErrorCode BSSCR_KSPLogConvergenceRate(PetscViewer v,PetscInt monitor_index, KSP ksp)
{
	PetscInt i;
	PetscReal *rlog;
	PetscLogDouble *log;
	PetscInt nits;
	PetscReal rho;

	BSSCR_KSPLogGetResidualTimeHistory(ksp,monitor_index, &nits,&rlog,&log);
	PetscViewerASCIIPrintf( v, "#-----------------------------------------------------------------------------------------------------------------------\n");
	PetscViewerASCIIPrintf( v, "# it         |r|              |r_0|               |r|/|r_0|          -log10(r/r0)        rho     |rn|/|r0|^{1/n}        \n" );
	PetscViewerASCIIPrintf( v, "#-----------------------------------------------------------------------------------------------------------------------\n");
	rho = 0.0;
	i = 0;
	PetscViewerASCIIPrintf( v, "%1.4d %14.12e %14.12e %14.12e %14.12e %14.12e %14.12e\n", i, rlog[i],rlog[0],rlog[i]/rlog[0],-log10(rlog[i]/rlog[0]), rho, 0.0 );
	for( i=1; i<nits; i++ ) {
		rho = (-log10(rlog[i]/rlog[0]))-(-log10(rlog[i-1]/rlog[0]));
		PetscViewerASCIIPrintf( v, "%1.4d %14.12e %14.12e %14.12e %14.12e %14.12e %14.12e \n",
				i, rlog[i],rlog[0],rlog[i]/rlog[0],-log10(rlog[i]/rlog[0]), rho, pow(rlog[i]/rlog[0], 1.0/(double)i) );
	}

	PetscFunctionReturn(0);
}


text:
#include <stdio.h>      //printf
#include <string.h>     //memset
#include <time.h>       //time
#include <cblas.h>      //cblas_scopy, cblas_sscal, cblas_sgemm
#include <sys/times.h>  //times
#include <unistd.h>     //sysconf
#include <math.h>       //fabsf
#include <xmmintrin.h>  //_mm_malloc
#define  ALIGN 16

float* gen_matrix(int order, float range){
    float* m = _mm_malloc(order * order * sizeof(float), ALIGN);
    srand(time(NULL));
    for(int i = 0; i < order * order; ++i)
        m[i] = ((float)rand() / (float)(RAND_MAX)) * range;
    return m;
}

void transpose_matrix(float* m, int order){
    for(int i = 0; i < order; ++i)
        for(int j = i; j < order; ++j){
            float tmp = m[i * order + j];
            m[i * order + j] = m[j * order + i];
            m[j * order + i] = tmp;
        }
}

void print_matrix(float* m, int order, FILE* out){
    fprintf(out, "Order: %d\n", order);
    for(int i = 0; i < order; ++i){
        for(int j = 0; j < order; ++j)
            fprintf(out, "%.4f ", m[i * order + j]);
        fprintf(out, "\n");
    }
    fprintf(out, "\n");
}

void get_matrix_norms(float* m, int order, float* l1_norm, float* max_norm){
    *l1_norm = 0;
    *max_norm = 0;
    for(int i = 0; i < order; ++i){
        float row_sum = 0;
        float col_sum = 0;
        for(int j = 0; j < order; ++j){
            row_sum += fabs(m[i * order + j]);
            col_sum += fabs(m[j * order + i]);
        }
        if(*l1_norm < col_sum)
            *l1_norm = col_sum;
        if(*max_norm < row_sum)
            *max_norm = row_sum;
    }
}

float* get_identity_matrix(int order){
    float* m = _mm_malloc(order * order * sizeof(float), ALIGN);
    memset(m, 0, order * order * sizeof(float));
    for(int i = 0; i < order; ++i)
        m[i * order + i] = 1;
    return m;
}

float* invert_matrix(float* A, int order, int iter_number){
    float* B = _mm_malloc(order * order * sizeof(float), ALIGN);
    cblas_scopy(order * order, A, 1, B, 1);
    transpose_matrix(B, order);
    float l1_norm, max_norm;
    get_matrix_norms(A, order, &l1_norm, &max_norm);
    cblas_sscal(order * order, 1.0 / (l1_norm * max_norm), B, 1);
    float* R = get_identity_matrix(order);
    cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, order, order, order, -1, B, order, A, order, 1, R, order);
    float* R_1_deg = _mm_malloc(order * order * sizeof(float), ALIGN);
    cblas_scopy(order * order, R, 1, R_1_deg, 1);
    float* tmp = NULL;
    float* inv_A = get_identity_matrix(order);
    for(int i = 0; i < iter_number; ++i){
        tmp = _mm_malloc(order * order * sizeof(float), ALIGN);
        cblas_saxpy(order * order, 1, R, 1, inv_A, 1);
        cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, order, order, order, 1, R, order, R_1_deg, order, 0, tmp,
                    order);
        _mm_free(R);
        R = tmp;
    }
    cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, order, order, order, 1, inv_A, order, B, order, 0.0, tmp,
            order);
    _mm_free(inv_A);
    inv_A = tmp;
    _mm_free(B);
    _mm_free(R_1_deg);
    return inv_A;
}

int main(){
    FILE* input = fopen("input.txt", "r");
    FILE* output = fopen("output.txt", "w");
    int N = 0, M = 0;
    fscanf(input, "%d%d", &N, &M);
    float* A = gen_matrix(N, 5.0);
    time_t start_real = time(NULL);
    float* inv_A = invert_matrix(A, N, M);
    time_t finish_real = time(NULL);
    float* rez = _mm_malloc(N * N * sizeof(float), ALIGN);
    cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, N, N, N, 1, inv_A, N, A, N, 0, rez, N);
    print_matrix(rez, N, output);
    printf("Total real time: %ld sec.\n", finish_real - start_real);
    _mm_free(A);
    _mm_free(inv_A);
    _mm_free(rez);
    fclose(input);
    fclose(output);
    return 0;
}

max_stars_repo_path: The 2nd course/Computers and devices/laboratory work #7/with_cblas/cblas.c


text:
/*==========================================================
 * arrayStxmWithMask.c - example in MATLAB External Interfaces
 *
 * Multiplies an input scalar (multiplier)
 * times a 1xN matrix (inMatrix)
 * and outputs a 1xN matrix (outMatrix)
 *
 * The calling syntax is:
 *
 *		outMatrix = arrayStxmWithMask(multiplier, inMatrix)
 *
 * This is a MEX-file for MATLAB.
 * Copyright 2007-2012 The MathWorks, Inc.
 *
 *========================================================*/

/* Compiled with
 * mex -lgsl -lgslcblas ./toolbox/nanodiffraction/mex/mex_pca.c
 */

#include "mex.h"
/* #include <stdio.h> */
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>

#define PI 3.14159265

/* The computational routine */
void pca(double *x, double *mask, double *qx, double *qy, double *w, double *phi, mwSize m, mwSize n)
{
    /*
    % first order moments
    tr = sum(dat(:));

    % calculate covariance
    cov_qx_qx= sum(sum(dat.*qy.^2))./tr;
    cov_mix  = sum(sum(dat.*qy.*qz))./tr; % Mischterm
    cov_qy_qy= sum(sum(dat.*qz.^2))./tr;

    % covariance matrix
    A=[cov_qx_qx cov_mix; cov_mix cov_qy_qy];

    % solve eigenvalue problem (V: Eigenvectoren, D:
    % Diagonalmatrix)
    [V,D] = eig(A);
    */

    mwSize col;
    mwSize row;

    double tr = 0;

    /* total sum of x */
    for (row=0; row<m; row++) {
        for (col=0; col<n; col++) {
            x[row*m + col] = x[row*m + col] * mask[row*m + col];
            tr += x[row*m + col];
        }
    }

    /* covariance matrix */
    double cov_qx = 0;
    double cov_mixed = 0;
    double cov_qy = 0;

    /* total sum of x */
    for (row=0; row<m; row++) {
        for (col=0; col<n; col++) {
            cov_qx += x[row*m + col] * pow(qx[row*m + col],2);
            cov_qy += x[row*m + col] * pow(qy[row*m + col],2);
            cov_mixed += x[row*m + col] * qx[row*m + col] * qy[row*m + col];
        }
    }

    double data[] = {cov_qx/tr   , cov_mixed/tr,
                     cov_mixed/tr, cov_qy/tr   };

    gsl_matrix_view mv
        = gsl_matrix_view_array (data, 2, 2);

    gsl_vector *eval = gsl_vector_alloc (2);
    gsl_matrix *evec = gsl_matrix_alloc (2, 2);

    gsl_eigen_symmv_workspace * ws
        = gsl_eigen_symmv_alloc (2);

    gsl_eigen_symmv (&mv.matrix, eval, evec, ws);

    gsl_eigen_symmv_free (ws);

    gsl_eigen_symmv_sort (eval, evec,
                        GSL_EIGEN_SORT_ABS_DESC);

    /* get eigenvalues and eigenvectors (assuming they are sorted) */
    double eval_i = gsl_vector_get (eval, 0);
    double eval_j = gsl_vector_get (eval, 1);

    gsl_vector_view evec_i
       = gsl_matrix_column (evec, 0);

    double xcomp = gsl_vector_get (&evec_i.vector, 0);
    double ycomp = gsl_vector_get (&evec_i.vector, 1);

    /* anisotropy */
    *w = (eval_i - eval_j) / (eval_i + eval_j);

    /* orientation */
    *phi = atan2(ycomp,xcomp) * 180.0 / PI;
    if (*phi < 0){
        *phi += 180;
    }

    gsl_vector_free (eval);
    gsl_matrix_free (evec);

    /* debugging information */
    /*
    printf ("%g\t%g\t%g\n", cov_qx/tr, cov_mixed/tr, cov_qy/tr );

    printf ("eigenvalue 1 = %g\n", eval_i);
    printf ("eigenvalue 2 = %g\n", eval_j);

    printf ("eigenvector = \n");
    gsl_vector_fprintf (stdout,
                        &evec_i.vector, "%g");

    printf ("x = %g\n", xcomp);
    printf ("y = %g\n", ycomp);
    printf ("w = %g\n", *w);
    printf ("phi = %g\n", *phi);
    */
}

/* The gateway function */
void mexFunction( int nlhs, mxArray *plhs[],        /* output variables */
                  int nrhs, const mxArray *prhs[])  /* input variables */
{
    /* double placeholder; */   /* input scalar */
    double *inData;         /* 1xN input matrix */
    double *inMask;         /* 1xN input matrix */
    double *qx;             /* 1xN input matrix */
    double *qy;             /* 1xN input matrix */
    size_t M;               /* size of matrix */
    size_t N;               /* size of matrix */
    double *w;              /* output scalar or matrix */
    double *phi;            /* output scalar or matrix */

    /* check for proper number of arguments */
    if(nrhs!=4) {
        mexErrMsgIdAndTxt("MyToolbox:arrayStxmWithMask:nrhs","Four inputs required.");
    }
    if(nlhs!=2) {
        mexErrMsgIdAndTxt("MyToolbox:arrayStxmWithMask:nlhs","Two outputs required.");
    }

    /* get the value of the scalar input  */
    /* qrMin = mxGetScalar(prhs[1]); */

    /* create a pointer to the data */
    inData  = mxGetPr(prhs[0]);
    inMask  = mxGetPr(prhs[1]);
    qx      = mxGetPr(prhs[2]);
    qy      = mxGetPr(prhs[3]);

    /* get dimensions of the input data matrix */
    M = mxGetM(prhs[0]);
    N = mxGetN(prhs[0]);

    /* create the output matrix */
    plhs[0] = mxCreateDoubleScalar(0);
    plhs[1] = mxCreateDoubleScalar(0);

    /* get a pointer to the real data in the output matrix */
    w   = mxGetPr(plhs[0]);
    phi = mxGetPr(plhs[1]);

    /* call the computational routine */
    pca(inData,inMask,qx,qy,w,phi,(mwSize)M,(mwSize)N);

}

max_stars_repo_name: irpgoe/nanodiffraction

max_issues_repo_name: irpgoe/nanodiffraction

max_forks_repo_name: irpgoe/nanodiffraction


text:
#ifndef NETWORK_CANUDPRECEIVER_H
#define NETWORK_CANUDPRECEIVER_H

#include <cstdint>
#include <gsl/gsl>
#include <QObject>

#include "tincan/canrawframe.h"
#include "network/udpasyncreceiver.h"

namespace network {

class Can_udp_receiver final : public QObject, public udp::Async_receiver
{
  Q_OBJECT

public:
  Can_udp_receiver() = default;
  Can_udp_receiver(const Can_udp_receiver&) = delete;
  Can_udp_receiver(Can_udp_receiver&&) = delete;
  Can_udp_receiver& operator=(const Can_udp_receiver&) = delete;
  Can_udp_receiver& operator=(Can_udp_receiver&&) = delete;

  void handle_receive(gsl::span<std::uint8_t> buffer) override;

signals:
  void received_frame(std::uint64_t, tin::Can_raw_frame);
};

}  // namespace network

#endif  // NETWORK_CANUDPRECEIVER_H

ext: h

max_stars_repo_name: mwkpe/tincan

max_issues_repo_name: jwkpeter/tincan

max_forks_repo_name: jwkpeter/tincan


text:
#include "utils.h"
#include <assert.h>
#include "stdio.h"
#include "string.h"
#include "cosmosis/datablock/c_datablock.h"
#include <gsl/gsl_math.h>
#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.h>

gsl_spline * spline_from_arrays(int n, double * x, double *y)
{
	gsl_spline * output = gsl_spline_alloc(gsl_interp_akima, n);
	assert (output!=NULL);
	gsl_spline_init(output,x,y,n);
	return output;
}

DATABLOCK_STATUS save_spline(c_datablock * block, const char * section,
	const char * n_name, const char * x_name, const char * y_name,
	gsl_spline * s)
{

	DATABLOCK_STATUS status = 0;
	if (strlen(n_name)>0){
		status |= c_datablock_put_int(block, section, n_name, s->size);
	}
	if (strlen(x_name)>0){
		status |= c_datablock_put_double_array_1d(block, section, x_name,
			s->x, s->size);
	}
	status |= c_datablock_put_double_array_1d(block, section, y_name,
		s->y, s->size);
	return status;
}

void reverse(double * x, int n)
{
	for (int i=0; i<n/2; i++){
		double tmp = x[i];
		x[i] = x[n-1-i];
		x[n-1-i] = tmp;
	}
}

static
double identity_function(double x, double y, double z, void * a)
{
  return z;
}

Interpolator2D *
load_interpolator_chi_function(c_datablock * block, gsl_spline * chi_of_z_spline,
	const char * section,
	const char * k_name, const char * z_name, const char * P_name,
	interp2d_modifier_function function, void * args
	)
{
	int nk, nz, nP;
	double *k=NULL, *z=NULL;
	double **P=NULL;
	int status = 0;

	status = c_datablock_get_double_grid(block, section,
		k_name, &nk, &k,
		z_name, &nz, &z,
		P_name, &P);

	if (status){
		fprintf(stderr, "Could not load interpolator for P(k).  Error %d\n",status);
		return NULL;
	}

        //read only P(k,z) and not P(k,0) since the growth factor will be inside the kernels
	for (int j=0; j<nk; j++){
		for (int i=0; i<nz; i++){
			P[j][i] = function(k[j], z[0], P[j][0], args);
		}
	}

	// What we have now is P(k, z).
	// We can optionally convert to P(k, chi)
	// If so we loop, lookup, and replace
	if (chi_of_z_spline){
		for (int i=0; i<nz; i++){
			double zi = z[i];
			double chi_i = gsl_spline_eval(chi_of_z_spline, zi, NULL);
			z[i] = chi_i;
		}
	}

	if (status) return NULL;
	Interpolator2D * interp = init_interp_2d_akima_grid(k, z, P, nk, nz);
	deallocate_2d_double(&P, nk);
	free(k);
	free(z);
	return interp;
}

Interpolator2D *
load_interpolator_chi(c_datablock * block, gsl_spline * chi_of_z_spline,
	const char * section,
	const char * k_name, const char * z_name, const char * P_name)
{

	return load_interpolator_chi_function(block, chi_of_z_spline, section, k_name, z_name, P_name, identity_function, NULL);
}

Interpolator2D *
load_interpolator_function(c_datablock * block,
	const char * section,
	const char * k_name, const char * z_name, const char * P_name,
	interp2d_modifier_function function, void * args
	)
{
	return load_interpolator_chi_function(block, NULL, section, k_name, z_name, P_name, function, args);

}

Interpolator2D *
load_interpolator(c_datablock * block,
	const char * section,
	const char * k_name, const char * z_name, const char * P_name)
{
	return load_interpolator_chi_function(block, NULL, section, k_name, z_name, P_name, identity_function, NULL);
}

// this is for the scale dependent growth rate from here----------------

Interpolator2D *
load_interpolator_chi_function2(c_datablock * block, gsl_spline * chi_of_z_spline,  gsl_spline * z_of_chi_spline,
	const char * section,
	const char * k_name, const char * z_name, const char * P_name,
	interp2d_modifier_function function, void * args
	)
{
	int nk, nz, nT;
	double *k=NULL, *z=NULL;
	double **P=NULL;
	int status = 0;

	status = c_datablock_get_double_grid(block, section,
		k_name, &nk, &k,
		z_name, &nz, &z,
		P_name, &P);

	double a[nk][nz],*pp[nk],**dTdz=NULL;
        for(int j=0;j<nk;j++){
	pp[j]=&a[j][0];}//for the derivative
	dTdz=pp;

	double b[nk][nz],*ppp[nk],**f=NULL;
        for(int j=0;j<nk;j++){
	ppp[j]=&b[j][0];}//for f(k,z)
	f=ppp;

	double c[nk][nz],*pppp[nk],**T=NULL;
        for(int j=0;j<nk;j++){
	pppp[j]=&c[j][0];}//for the T(k,z)
	T=pppp;

	if (status){
		fprintf(stderr, "Could not load interpolator for T(k).  Error %d\n",status);
		return NULL;
	}

	for (int j=0; j<nk; j++){
		for (int i=0; i<nz; i++){
			P[j][i] = function(k[j], z[i], P[j][i], args);
		}
	}

	//this is to take T(k,z)
	for (int j=0; j<nk; j++){
		for (int i=0; i<nz; i++){
			T[j][i] = sqrt(P[j][i]/P[j][0]);
		}
	}

        //derivative dT/dz

	for (int j=0; j<nk; j++){
		for (int i=0; i<nz; i++){

   	        //forward differences

   		if(i==0){
   			dTdz[j][i] = (T[j][i+1]-T[j][i])/(z[i+1]-z[i]);
			}

		//central differences

		if(0<i && i<(nz-1)){
   			dTdz[j][i] = (T[j][i+1]-T[j][i-1])/(2.*(z[i+1]-z[i]));
			}
		//backward differences

		if(i==(nz-1)){
   			dTdz[j][i] = (T[j][i]-T[j][i-1])/(z[i]-z[i-1]);
			}

		}
	}

	//printf("ok with derivative\n");

	//this is to take f(k,z)
	for (int j=0; j<nk; j++){
		for (int i=0; i<nz; i++){
			f[j][i] = -((z[i]+1.)/T[j][i])*dTdz[j][i];
		}
	}

	// What we have now is f(k, z).
	// We can optionally convert to f(k, z(chi))
	// If so we loop, lookup, and replace
	if (chi_of_z_spline){
		for (int i=0; i<nz; i++){
                        //int j=0; //my modification
			double zi = z[i]; //// j instead of i
			double chi_i = gsl_spline_eval(chi_of_z_spline, zi, NULL);
			//double z_i = gsl_spline_eval(z_of_chi_spline, chi_i, NULL);
			z[i] = chi_i;
		}
	}

	Interpolator2D * interp2 = init_interp_2d_akima_grid(k, z, f, nk, nz);
	return interp2;
}

Interpolator2D *
load_interpolator_chi2(c_datablock * block, gsl_spline * chi_of_z_spline, gsl_spline * z_of_chi_spline,
	const char * section,
	const char * k_name, const char * z_name, const char * P_name)
{

	return load_interpolator_chi_function2(block, chi_of_z_spline, z_of_chi_spline, section, k_name, z_name, P_name, identity_function, NULL);
}

Interpolator2D *
load_interpolator_function2(c_datablock * block,
	const char * section,
	const char * k_name, const char * z_name, const char * P_name,
	interp2d_modifier_function function, void * args
	)
{
	return load_interpolator_chi_function2(block, NULL, NULL, section, k_name, z_name, P_name, function, args);

}

Interpolator2D *
load_interpolator2(c_datablock * block,
	const char * section,
	const char * k_name, const char * z_name, const char * P_name)
{
	return load_interpolator_chi_function2(block, NULL, NULL, section, k_name, z_name, P_name, identity_function, NULL);
}

// to here---------------------------------------

// here write the scale dependent growth factor

Interpolator2D *
load_interpolator_chi_function3(c_datablock * block, gsl_spline * chi_of_z_spline,  gsl_spline * z_of_chi_spline,
	const char * section,
	const char * k_name, const char * z_name, const char * P_name,
	interp2d_modifier_function function, void * args
	)
{
	int nk, nz, nT;
	double *k=NULL, *z=NULL;
	double **P=NULL;
	int status = 0;

	status = c_datablock_get_double_grid(block, section,
		k_name, &nk, &k,
		z_name, &nz, &z,
		P_name, &P);

	double b[nk][nz],*ppp[nk],**D=NULL;
        for(int j=0;j<nk;j++){
	ppp[j]=&b[j][0];}//for D(k,z)
	D=ppp;

	if (status){
		fprintf(stderr, "Could not load interpolator for P_nl(k).  Error %d\n",status);
		return NULL;
	}

	for (int j=0; j<nk; j++){
		for (int i=0; i<nz; i++){
			P[j][i] = function(k[j], z[i], P[j][i], args);
		}
	}

	//printf("ok with derivative\n");

	//this is to take D(k,z)
	for (int j=0; j<nk; j++){
		for (int i=0; i<nz; i++){
			D[j][i] = sqrt(P[j][i]/P[j][0]);
		}
	}

	// What we have now is D(k, z).
	// We can optionally convert to D(k, z(chi))
	// If so we loop, lookup, and replace
	if (chi_of_z_spline){
		for (int i=0; i<nz; i++){
                        //int j=0; //my modification
			double zi = z[i]; //// j instead of i
			double chi_i = gsl_spline_eval(chi_of_z_spline, zi, NULL);
			//double z_i = gsl_spline_eval(z_of_chi_spline, chi_i, NULL);
			z[i] = chi_i;
		}
	}

	//printf("ok with z_of_chi\n");
	Interpolator2D * interp3 = init_interp_2d_akima_grid(k, z, D, nk, nz);
	return interp3;
}

Interpolator2D *
load_interpolator_chi3(c_datablock * block, gsl_spline * chi_of_z_spline, gsl_spline * z_of_chi_spline,
	const char * section,
	const char * k_name, const char * z_name, const char * P_name)
{

	return load_interpolator_chi_function3(block, chi_of_z_spline, z_of_chi_spline, section, k_name, z_name, P_name, identity_function, NULL);
}

Interpolator2D *
load_interpolator_function3(c_datablock * block,
	const char * section,
	const char * k_name, const char * z_name, const char * P_name,
	interp2d_modifier_function function, void * args
	)
{
	return load_interpolator_chi_function3(block, NULL, NULL, section, k_name, z_name, P_name, function, args);

}

Interpolator2D *
load_interpolator3(c_datablock * block,
	const char * section,
	const char * k_name, const char * z_name, const char * P_name)
{
	return load_interpolator_chi_function3(block, NULL, NULL, section, k_name, z_name, P_name, identity_function, NULL);
}

//-------------------to here

//-----------------here is the scale dependent galaxy bias accounting for neutrinos

Interpolator2D *
load_interpolator_chi_function4(c_datablock * block, gsl_spline * chi_of_z_spline,  gsl_spline * z_of_chi_spline,
	const char * section,
	const char * k_name, const char * z_name, const char * P_name,
	interp2d_modifier_function function, void * args
	)
{
	int nk, nz, nT;
	double *k=NULL, *z=NULL;
	double **P=NULL;
	int status = 0;

	status = c_datablock_get_double_grid(block, section,
		k_name, &nk, &k,
		z_name, &nz, &z,
		P_name, &P);

	double bb[nk][nz],*pppp[nk],**Bias=NULL;
        for(int j=0;j<nk;j++){
	pppp[j]=&bb[j][0];}//for b(k,z)
	Bias=pppp;

	double bbL[nk][nz],*ppppL[nk],**bL=NULL;
        for(int j=0;j<nk;j++){
	ppppL[j]=&bbL[j][0];}//for bL(k,z)
	bL=ppppL;

	double bbS[nk][nz],*ppppS[nk],**bS=NULL;
        for(int j=0;j<nk;j++){
	ppppS[j]=&bbS[j][0];}//for bS(k,z)
	bS=ppppS;

	if (status){
		fprintf(stderr, "Could not load interpolator for P_nl(k).  Error %d\n",status);
		return NULL;
	}

	for (int j=0; j<nk; j++){
		for (int i=0; i<nz; i++){
			P[j][i] = function(k[j], z[i], P[j][i], args);
		}
	}

  	int status7=0;
	double omega_m;
	status7 = c_datablock_get_double(block, COSMOLOGICAL_PARAMETERS_SECTION, "omega_m", &omega_m);
	double Om = omega_m;

  	int status8=0;
	double omega_vh2;
	status8 = c_datablock_get_double(block, COSMOLOGICAL_PARAMETERS_SECTION, "omnuh2", &omega_vh2);
	double Ovh2 = omega_vh2;

  	int status9=0;
	double hs;
	status9 = c_datablock_get_double(block, COSMOLOGICAL_PARAMETERS_SECTION, "h0", &hs);
	double hubble = hs;

        double fv = (Ovh2/(hubble * hubble))/Om;

        double mv = Ovh2 * 93.14; //this is Sum(mv)/1eV in h/Mpc but for one massive neutrino only Sum(mv)=mv

        // this the neutrino free streaming scale in h/Mpc
	//for (int i=0; i<nz; i++){
	//	kfs[0][i] = 0.82 * sqrt(1.0 - Om + Om * (1.0 + z[i]) * (1.0 + z[i]) * (1.0 + z[i]))/((1.0 + z[i]) * (1.0 + z[i])) * mv;
	//}
        double knr;
        knr = 0.018 * sqrt(Om) * sqrt(mv);

        double constant=1.0;
	//this is to take b(k,z)
	for (int j=0; j<nk; j++){
		for (int i=0; i<nz; i++){

                bL[j][i] =  constant/sqrt(P[0][i]/P[0][0]);

                bS[j][i] = constant/sqrt(P[0][i]/P[0][0]) - constant/sqrt(P[0][i]/P[0][0]) * fv;

                Bias[j][i] = bL[j][i] + (bS[j][i] - bL[j][i]) * 0.5 * (tanh(log(pow(k[j]/(knr),5.0)))+1.0) ;

		}
	}

	// What we have now is b(k, z).
	// We can optionally convert to b(k, z(chi))
	// If so we loop, lookup, and replace
	if (chi_of_z_spline){
		for (int i=0; i<nz; i++){
                        //int j=0; //my modification
			double zi = z[i]; //// j instead of i
			double chi_i = gsl_spline_eval(chi_of_z_spline, zi, NULL);
			//double z_i = gsl_spline_eval(z_of_chi_spline, chi_i, NULL);
			z[i] = chi_i;
		}
	}

	//printf("ok with z_of_chi\n");
	Interpolator2D * interp4 = init_interp_2d_akima_grid(k, z, Bias, nk, nz);
	return interp4;
}

Interpolator2D *
load_interpolator_chi4(c_datablock * block, gsl_spline * chi_of_z_spline, gsl_spline * z_of_chi_spline,
	const char * section,
	const char * k_name, const char * z_name, const char * P_name)
{

	return load_interpolator_chi_function4(block, chi_of_z_spline, z_of_chi_spline, section, k_name, z_name, P_name, identity_function, NULL);
}

Interpolator2D *
load_interpolator_function4(c_datablock * block,
	const char * section,
	const char * k_name, const char * z_name, const char * P_name,
	interp2d_modifier_function function, void * args
	)
{
	return load_interpolator_chi_function4(block, NULL, NULL, section, k_name, z_name, P_name, function, args);

}

Interpolator2D *
load_interpolator4(c_datablock * block,
	const char * section,
	const char * k_name, const char * z_name, const char * P_name)
{
	return load_interpolator_chi_function4(block, NULL, NULL, section, k_name, z_name, P_name, identity_function, NULL);
}

//------------------to here

gsl_spline * load_spline(c_datablock * block, const char * section,
	const char * x_name, const char * y_name)
{
	int status = 0 ;
	int nx, ny;
	double *x, *y;
	status |= c_datablock_get_double_array_1d(block, section, x_name, &x, &nx);
	if (status) return NULL;
	status |= c_datablock_get_double_array_1d(block, section, y_name, &y, &ny);
	if (status){
		free(x);
		return NULL;
	}

	gsl_spline * output = NULL;

	if (x[1]<x[0]) {reverse(x, nx); reverse(y, ny);}

	if (nx==ny) output=spline_from_arrays(nx, x, y);
	else {fprintf(stderr, "Non-matching array sizes %s=%d, %s=%d\n", x_name, nx, y_name, ny);}
	free(x);
	free(y);
	return output;
}

/*
gsl_spline2d *
load_gsl_interpolator_chi(c_datablock * block, gsl_spline * chi_of_z_spline,
	const char * section,
	const char * k_name, const char * z_name, const char * P_name,
	)
{
	int nk, nz, nP;
	double *k=NULL, *z=NULL;
	double **P=NULL;
	int status = 0;

	status = c_datablock_get_double_grid(block, section,
		k_name, &nk, &k,
		z_name, &nz, &z,
		P_name, &P);

	if (status){
		fprintf(stderr, "Could not load interpolator for P(k).  Error %d\n",status);
		return NULL;
	}

	// What we have now is P(k, z).
	// We can optionally convert to P(k, chi)
	// If so we loop, lookup, and replace
	if (chi_of_z_spline){
		for (int i=0; i<nz; i++){
			double zi = z[i];
			double chi_i = gsl_spline_eval(chi_of_z_spline, zi, NULL);
			z[i] = chi_i;
		}
	}

	if (status) return NULL;
        const gsl_interp2d_type *T = gsl_interp2d_bicubic;
	gsl_spline2d *spline = gsl_spline2d_alloc(T, nz, nk);
	gsl_interp_accel *xacc = gsl_interp_accel_alloc();
	gsl_interp_accel *yacc = gsl_interp_accel_alloc();
	gsl_spline_init

	Interpolator2D * interp = init_interp_2d_akima_grid(k, z, P, nk, nz);
	deallocate_2d_double(&P, nk);
	free(k);
	free(z);
	return interp;
}
*/


text:
static char help[] =
"Unstructured 2D FEM solution of nonlinear Poisson problem.  Option prefix un_.\n"
"Equation is\n"
"    - div( a(u,x,y) grad u ) = f(u,x,y)\n"
"on arbitrary 2D polygonal domain, with boundary data g_D(x,y), g_N(x,y).\n"
"Functions a(), f(), g_D(), g_N(), and u_exact() are given as formulas.\n"
"Input files in PETSc binary format contain node coordinates, elements, Neumann\n"
"boundary segments and boundary flags.  Allows non-homogeneous Dirichlet and/or\n"
"Neumann boundary conditions.  Five different solution cases are implemented.\n\n";

#include <petsc.h>
#include "../quadrature.h"
#include "um.h"
#include "cases.h"

//STARTCTX
typedef struct {
    UM     *mesh;
    int    solncase,
           quaddegree;
    double (*a_fcn)(double, double, double);
    double (*f_fcn)(double, double, double);
    double (*gD_fcn)(double, double);
    double (*gN_fcn)(double, double);
    double (*uexact_fcn)(double, double);
    PetscLogStage readstage, setupstage, solverstage, resstage, jacstage;  //STRIP
} unfemCtx;
//ENDCTX

//STARTFEM
double chi(int L, double xi, double eta) {
    const double z[3] = {1.0 - xi - eta, xi, eta};
    return z[L];
}

const double dchi[3][2] = {{-1.0,-1.0},{ 1.0, 0.0},{ 0.0, 1.0}};

// evaluate v(xi,eta) on reference element using local node numbering
double eval(const double v[3], double xi, double eta) {
    double sum = 0.0;
    int    L;
    for (L = 0; L < 3; L++)
        sum += v[L] * chi(L,xi,eta);
    return sum;
}

double InnerProd(const double V[2], const double W[2]) {
    return V[0] * W[0] + V[1] * W[1];
}
//ENDFEM

extern PetscErrorCode FillExact(Vec, unfemCtx*);
extern PetscErrorCode FormFunction(SNES, Vec, Vec, void*);
extern PetscErrorCode FormPicard(SNES, Vec, Mat, Mat, void*);
extern PetscErrorCode Preallocation(Mat, unfemCtx*);

int main(int argc,char **argv) {
    PetscErrorCode ierr;
    PetscBool   viewmesh = PETSC_FALSE,
                viewsoln = PETSC_FALSE,
                noprealloc = PETSC_FALSE,
                savepintbinary = PETSC_FALSE,
                savepintmatlab = PETSC_FALSE;
    char        root[256] = "", nodesname[256], issname[256], solnname[256],
                pintname[256] = "";
    int         savepintlevel = -1, levels;
    UM          mesh;
    unfemCtx    user;
    SNES        snes;
    KSP         ksp;
    PC          pc;
    PCType      pctype;
    Mat         A;
    Vec         r, u, uexact;
    double      err, h_max;

    PetscInitialize(&argc,&argv,NULL,help);
    ierr = PetscLogStageRegister("Read mesh      ", &user.readstage); CHKERRQ(ierr);  //STRIP
    ierr = PetscLogStageRegister("Set-up         ", &user.setupstage); CHKERRQ(ierr);  //STRIP
    ierr = PetscLogStageRegister("Solver         ", &user.solverstage); CHKERRQ(ierr);  //STRIP
    ierr = PetscLogStageRegister("Residual eval  ", &user.resstage); CHKERRQ(ierr);  //STRIP
    ierr = PetscLogStageRegister("Jacobian eval  ", &user.jacstage); CHKERRQ(ierr);  //STRIP

    user.quaddegree = 1;
    user.solncase = 0;
    ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "un_", "options for unfem", ""); CHKERRQ(ierr);
    ierr = PetscOptionsInt("-case",
           "exact solution cases: 0=linear, 1=nonlinear, 2=nonhomoNeumann, 3=chapter3, 4=koch",
           "unfem.c",user.solncase,&(user.solncase),NULL); CHKERRQ(ierr);
    ierr = PetscOptionsString("-gamg_save_pint_binary",
           "filename under which to save interpolation operator (Mat) in PETSc binary format",
           "unfem.c",pintname,pintname,sizeof(pintname),&savepintbinary); CHKERRQ(ierr);
    ierr = PetscOptionsString("-gamg_save_pint_matlab",
           "filename under which to save interpolation operator (Mat) in ascii Matlab format",
           "unfem.c",pintname,pintname,sizeof(pintname),&savepintmatlab); CHKERRQ(ierr);
    ierr = PetscOptionsInt("-gamg_save_pint_level",
           "saved interpolation operator is between L-1 and L where this option sets L; defaults to finest levels",
           "unfem.c",savepintlevel,&savepintlevel,NULL); CHKERRQ(ierr);
    ierr = PetscOptionsString("-mesh",
           "file name root of mesh stored in PETSc binary with .vec,.is extensions",
           "unfem.c",root,root,sizeof(root),NULL); CHKERRQ(ierr);
    ierr = PetscOptionsBool("-noprealloc",
           "do not perform preallocation before matrix assembly",
           "unfem.c",noprealloc,&noprealloc,NULL); CHKERRQ(ierr);
    ierr = PetscOptionsInt("-quaddegree",
           "quadrature degree (1,2,3)",
           "unfem.c",user.quaddegree,&(user.quaddegree),NULL); CHKERRQ(ierr);
    ierr = PetscOptionsBool("-view_mesh",
           "view loaded mesh (nodes and elements) at stdout",
           "unfem.c",viewmesh,&viewmesh,NULL); CHKERRQ(ierr);
    ierr = PetscOptionsBool("-view_solution",
           "view solution u(x,y) to binary file; uses root name of mesh plus .soln\nsee petsc2tricontour.py to view graphically",
           "unfem.c",viewsoln,&viewsoln,NULL); CHKERRQ(ierr);
    ierr = PetscOptionsEnd(); CHKERRQ(ierr);

    // determine filenames
    if (strlen(root) == 0) {
        SETERRQ(PETSC_COMM_WORLD,4,"no mesh name root given; rerun with '-un_mesh foo'");
    }
    strcpy(nodesname, root);
    strncat(nodesname, ".vec", 4);
    strcpy(issname, root);
    strncat(issname, ".is", 3);

    // set source/boundary functions and exact solution
    user.a_fcn = &a_lin;
    user.f_fcn = &f_lin;
    user.uexact_fcn = &uexact_lin;
    user.gD_fcn = &gD_lin;
    user.gN_fcn = &gN_lin;
    switch (user.solncase) {
        case 0 :
            break;
        case 1 :
            user.a_fcn = &a_nonlin;
            user.f_fcn = &f_nonlin;
            break;
        case 2 :
            user.gN_fcn = &gN_linneu;
            break;
        case 3 :
            user.a_fcn = &a_square;
            user.f_fcn = &f_square;
            user.uexact_fcn = &uexact_square;
            user.gD_fcn = &gD_square;
            user.gN_fcn = NULL;  // seg fault if ever called
            break;
        case 4 :
            user.a_fcn = &a_koch;
            user.f_fcn = &f_koch;
            user.uexact_fcn = NULL;  // seg fault if ever called
            user.gD_fcn = &gD_koch;
            user.gN_fcn = NULL;  // seg fault if ever called
            break;
        default :
            SETERRQ(PETSC_COMM_WORLD,1,"other solution cases not implemented");
    }

    PetscLogStagePush(user.readstage);
//STARTREADMESH
    // read mesh object of type UM
    ierr = UMInitialize(&mesh); CHKERRQ(ierr);
    ierr = UMReadNodes(&mesh,nodesname); CHKERRQ(ierr);
    ierr = UMReadISs(&mesh,issname); CHKERRQ(ierr);
    ierr = UMStats(&mesh, &h_max, NULL, NULL, NULL); CHKERRQ(ierr);
    user.mesh = &mesh;
//ENDREADMESH
    PetscLogStagePop();

    if (viewmesh) {
        PetscViewer stdoutviewer;
        ierr = PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&stdoutviewer); CHKERRQ(ierr);
        ierr = UMViewASCII(&mesh,stdoutviewer); CHKERRQ(ierr);
    }

    PetscLogStagePush(user.setupstage);
//STARTMAININITIAL
    // configure Vecs and SNES
    ierr = VecCreate(PETSC_COMM_WORLD,&r); CHKERRQ(ierr);
    ierr = VecSetSizes(r,PETSC_DECIDE,mesh.N); CHKERRQ(ierr);
    ierr = VecSetFromOptions(r); CHKERRQ(ierr);
    ierr = VecDuplicate(r,&u); CHKERRQ(ierr);
    ierr = VecSet(u,0.0); CHKERRQ(ierr);
    ierr = SNESCreate(PETSC_COMM_WORLD,&snes); CHKERRQ(ierr);
    ierr = SNESSetFunction(snes,r,FormFunction,&user); CHKERRQ(ierr);

    // reset default KSP and PC
    ierr = SNESGetKSP(snes,&ksp); CHKERRQ(ierr);
    ierr = KSPSetType(ksp,KSPCG); CHKERRQ(ierr);
    ierr = KSPGetPC(ksp,&pc); CHKERRQ(ierr);
    ierr = PCSetType(pc,PCICC); CHKERRQ(ierr);

    // setup matrix for Picard iteration, including preallocation
    ierr = MatCreate(PETSC_COMM_WORLD,&A); CHKERRQ(ierr);
    ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,mesh.N,mesh.N); CHKERRQ(ierr);
    ierr = MatSetFromOptions(A); CHKERRQ(ierr);
    ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE); CHKERRQ(ierr);
    if (noprealloc) {
        ierr = MatSetUp(A); CHKERRQ(ierr);
    } else {
        ierr = Preallocation(A,&user); CHKERRQ(ierr);
    }
    ierr = SNESSetJacobian(snes,A,A,FormPicard,&user); CHKERRQ(ierr);
    ierr = SNESSetFromOptions(snes); CHKERRQ(ierr);
    PetscLogStagePop();  //STRIP

    // solve
    PetscLogStagePush(user.solverstage);  //STRIP
    ierr = SNESSolve(snes,NULL,u);CHKERRQ(ierr);
//ENDMAININITIAL
    PetscLogStagePop();

    // report if PC is GAMG
    ierr = PCGetType(pc,&pctype); CHKERRQ(ierr);
    if (strcmp(pctype,"gamg") == 0) {
        ierr = PCMGGetLevels(pc,&levels); CHKERRQ(ierr);
        ierr = PetscPrintf(PETSC_COMM_WORLD,
               "  PC is GAMG with %d levels\n",levels); CHKERRQ(ierr);
    }

    // save Pint from GAMG if requested
    if (strlen(pintname) > 0) {
        Mat         pint;
        PetscViewer viewer;
        if (strcmp(pctype,"gamg") != 0) {
            SETERRQ(PETSC_COMM_WORLD,2,"option -un_gamg_save_pint set but PC is not of type PCGAMG");
        }
        if (savepintlevel >= levels) {
            SETERRQ(PETSC_COMM_WORLD,3,"invalid level given in -un_gamg_save_pint_level");
        }
        if (savepintbinary && savepintmatlab) {
            SETERRQ(PETSC_COMM_WORLD,4,"only one of -un_gamg_save_pint_binary OR -un_gamg_save_pint_matlab is allowed");
        }
        if (savepintlevel <= 0) {
            savepintlevel = levels - 1;
        }
        ierr = PCMGGetInterpolation(pc,savepintlevel,&pint); CHKERRQ(ierr);
        if (savepintbinary) {
            ierr = PetscPrintf(PETSC_COMM_WORLD,
               "  saving interpolation operator on levels %d->%d in binary format to %s ...\n",
               savepintlevel-1,savepintlevel,pintname); CHKERRQ(ierr);
            ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,pintname,FILE_MODE_WRITE,&viewer); CHKERRQ(ierr);
            ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB); CHKERRQ(ierr);
        } else {
            ierr = PetscPrintf(PETSC_COMM_WORLD,
               "  saving interpolation operator on levels %d->%d in ascii format to %s ...\n",
               savepintlevel-1,savepintlevel,pintname); CHKERRQ(ierr);
            ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD,pintname,&viewer); CHKERRQ(ierr);
            ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB); CHKERRQ(ierr);
        }
        ierr = MatView(pint,viewer); CHKERRQ(ierr);
    }

    // if exact solution available, report numerical error
    if (user.uexact_fcn) {
        ierr = VecDuplicate(r,&uexact); CHKERRQ(ierr);
        ierr = FillExact(uexact,&user); CHKERRQ(ierr);
        ierr = VecAXPY(u,-1.0,uexact); CHKERRQ(ierr);    // u <- u + (-1.0) uexact
        ierr = VecNorm(u,NORM_INFINITY,&err); CHKERRQ(ierr);
        ierr = PetscPrintf(PETSC_COMM_WORLD,
                   "case %d result for N=%d nodes with h = %.3e :  |u-u_ex|_inf = %.2e\n",
                   user.solncase,mesh.N,h_max,err); CHKERRQ(ierr);
        VecDestroy(&uexact);
    } else {
        ierr = PetscPrintf(PETSC_COMM_WORLD,
                   "case %d result for N=%d nodes with h = %.3e ... done\n",
                   user.solncase,mesh.N,h_max); CHKERRQ(ierr);
    }

    // save solution in PETSc binary if requested
    if (viewsoln) {
        strcpy(solnname, root);
        strncat(solnname, ".soln", 5);
        ierr = PetscPrintf(PETSC_COMM_WORLD,
                   "writing solution in binary format to %s ...\n",solnname); CHKERRQ(ierr);
        ierr = UMViewSolutionBinary(&mesh,solnname,u); CHKERRQ(ierr);
    }

    // clean-up
    VecDestroy(&u);  VecDestroy(&r);
    MatDestroy(&A);  SNESDestroy(&snes);  UMDestroy(&mesh);
    return PetscFinalize();
}

PetscErrorCode FillExact(Vec uexact, unfemCtx *ctx) {
    PetscErrorCode ierr;
    const Node   *aloc;
    double       *auexact;
    int          i;
    ierr = UMGetNodeCoordArrayRead(ctx->mesh,&aloc); CHKERRQ(ierr);
    ierr = VecGetArray(uexact,&auexact); CHKERRQ(ierr);
    for (i = 0; i < ctx->mesh->N; i++) {
        auexact[i] = ctx->uexact_fcn(aloc[i].x,aloc[i].y);
    }
    ierr = VecRestoreArray(uexact,&auexact); CHKERRQ(ierr);
    ierr = UMRestoreNodeCoordArrayRead(ctx->mesh,&aloc); CHKERRQ(ierr);
    return 0;
}

//STARTRESIDUAL
PetscErrorCode FormFunction(SNES snes, Vec u, Vec F, void *ctx) {
    PetscErrorCode ierr;
    unfemCtx        *user = (unfemCtx*)ctx;
    const Quad2DTri q = symmgauss[user->quaddegree-1];
    const int       *ae, *ans, *abf, *en;
    const Node      *aloc;
    const double    *au;
    int             p, na, nb, k, l, r;
    double          *aF, unode[3], gradu[2], gradpsi[3][2], uquad[4],
                    aquad[4], fquad[4], dx, dy, dx1, dx2, dy1, dy2,
                    detJ, ls, xmid, ymid, sint, xx, yy, psi, ip, sum;

    PetscLogStagePush(user->resstage);  //STRIP
    ierr = VecSet(F,0.0); CHKERRQ(ierr);
    ierr = VecGetArray(F,&aF); CHKERRQ(ierr);
    ierr = UMGetNodeCoordArrayRead(user->mesh,&aloc); CHKERRQ(ierr);
    ierr = ISGetIndices(user->mesh->bf,&abf); CHKERRQ(ierr);

    // Neumann boundary segment contributions (if any)
    if (user->mesh->P > 0) {
        ierr = ISGetIndices(user->mesh->ns,&ans); CHKERRQ(ierr);
        for (p = 0; p < user->mesh->P; p++) {
            na = ans[2*p+0];  nb = ans[2*p+1];  // nodes at end of segment
            dx = aloc[na].x-aloc[nb].x;  dy = aloc[na].y-aloc[nb].y;
            ls = sqrt(dx * dx + dy * dy);  // length of segment
            // midpoint rule; psi_na=psi_nb=0.5 at midpoint of segment
            xmid = 0.5*(aloc[na].x+aloc[nb].x);
            ymid = 0.5*(aloc[na].y+aloc[nb].y);
            sint = 0.5 * ls * user->gN_fcn(xmid,ymid);
            // nodes could be Dirichlet
            if (abf[na] != 2)
                aF[na] -= sint;
            if (abf[nb] != 2)
                aF[nb] -= sint;
        }
        ierr = ISRestoreIndices(user->mesh->ns,&ans); CHKERRQ(ierr);
    }

    // element contributions and Dirichlet node residuals
    ierr = VecGetArrayRead(u,&au); CHKERRQ(ierr);
    ierr = ISGetIndices(user->mesh->e,&ae); CHKERRQ(ierr);
    for (k = 0; k < user->mesh->K; k++) {
        // element geometry and hat function gradients
        en = ae + 3*k;  // en[0], en[1], en[2] are nodes of element k
        dx1 = aloc[en[1]].x - aloc[en[0]].x;
        dx2 = aloc[en[2]].x - aloc[en[0]].x;
        dy1 = aloc[en[1]].y - aloc[en[0]].y;
        dy2 = aloc[en[2]].y - aloc[en[0]].y;
        detJ = dx1 * dy2 - dx2 * dy1;
        for (l = 0; l < 3; l++) {
            gradpsi[l][0] = ( dy2 * dchi[l][0] - dy1 * dchi[l][1]) / detJ;
            gradpsi[l][1] = (-dx2 * dchi[l][0] + dx1 * dchi[l][1]) / detJ;
        }
        // u and grad u on element
        gradu[0] = 0.0;
        gradu[1] = 0.0;
        for (l = 0; l < 3; l++) {
            if (abf[en[l]] == 2)  // enforces symmetry
                unode[l] = user->gD_fcn(aloc[en[l]].x,aloc[en[l]].y);
            else
                unode[l] = au[en[l]];
            gradu[0] += unode[l] * gradpsi[l][0];
            gradu[1] += unode[l] * gradpsi[l][1];
        }
        // function values at quadrature points on element
        for (r = 0; r < q.n; r++) {
            uquad[r] = eval(unode,q.xi[r],q.eta[r]);
            xx = aloc[en[0]].x + dx1 * q.xi[r] + dx2 * q.eta[r];
            yy = aloc[en[0]].y + dy1 * q.xi[r] + dy2 * q.eta[r];
            aquad[r] = user->a_fcn(uquad[r],xx,yy);
            fquad[r] = user->f_fcn(uquad[r],xx,yy);
        }
        // residual contribution for each node of element
        for (l = 0; l < 3; l++) {
            if (abf[en[l]] == 2) { // set Dirichlet residual
                xx = aloc[en[l]].x;   yy = aloc[en[l]].y;
                aF[en[l]] = au[en[l]] - user->gD_fcn(xx,yy);
            } else {
                sum = 0.0;
                for (r = 0; r < q.n; r++) {
                    psi = chi(l,q.xi[r],q.eta[r]);
                    ip  = InnerProd(gradu,gradpsi[l]);
                    sum += q.w[r] * ( aquad[r] * ip - fquad[r] * psi );
                }
                aF[en[l]] += fabs(detJ) * sum;
            }
        }
    }

    ierr = ISRestoreIndices(user->mesh->e,&ae); CHKERRQ(ierr);
    ierr = VecRestoreArrayRead(u,&au); CHKERRQ(ierr);
    ierr = ISRestoreIndices(user->mesh->bf,&abf); CHKERRQ(ierr);
    ierr = UMRestoreNodeCoordArrayRead(user->mesh,&aloc); CHKERRQ(ierr);
    ierr = VecRestoreArray(F,&aF); CHKERRQ(ierr);
    PetscLogStagePop();  //STRIP
    return 0;
}
//ENDRESIDUAL

PetscErrorCode FormPicard(SNES snes, Vec u, Mat A, Mat P, void *ctx) {
    PetscErrorCode ierr;
    unfemCtx        *user = (unfemCtx*)ctx;
    const Quad2DTri q = symmgauss[user->quaddegree-1];
    const int       *ae, *abf, *en;
    const Node      *aloc;
    const double    *au;
    double          unode[3], gradpsi[3][2], uquad[4], aquad[4], v[9],
                    dx1, dx2, dy1, dy2, detJ, xx, yy, sum;
    int             n, k, l, m, r, cr, cv, row[3];

    PetscLogStagePush(user->jacstage);  //STRIP
    ierr = MatZeroEntries(P); CHKERRQ(ierr);
    ierr = ISGetIndices(user->mesh->bf,&abf); CHKERRQ(ierr);
    for (n = 0; n < user->mesh->N; n++) {
        if (abf[n] == 2) {
            v[0] = 1.0;
            ierr = MatSetValues(P,1,&n,1,&n,v,ADD_VALUES); CHKERRQ(ierr);
        }
    }
    ierr = ISGetIndices(user->mesh->e,&ae); CHKERRQ(ierr);
    ierr = VecGetArrayRead(u,&au); CHKERRQ(ierr);
    ierr = UMGetNodeCoordArrayRead(user->mesh,&aloc); CHKERRQ(ierr);
    for (k = 0; k < user->mesh->K; k++) {
        en = ae + 3*k;  // en[0], en[1], en[2] are nodes of element k
        // geometry of element
        dx1 = aloc[en[1]].x - aloc[en[0]].x;
        dx2 = aloc[en[2]].x - aloc[en[0]].x;
        dy1 = aloc[en[1]].y - aloc[en[0]].y;
        dy2 = aloc[en[2]].y - aloc[en[0]].y;
        detJ = dx1 * dy2 - dx2 * dy1;
        // gradients of hat functions and u on element
        for (l = 0; l < 3; l++) {
            gradpsi[l][0] = ( dy2 * dchi[l][0] - dy1 * dchi[l][1]) / detJ;
            gradpsi[l][1] = (-dx2 * dchi[l][0] + dx1 * dchi[l][1]) / detJ;
            if (abf[en[l]] == 2)
                unode[l] = user->gD_fcn(aloc[en[l]].x,aloc[en[l]].y);
            else
                unode[l] = au[en[l]];
        }
        // function values at quadrature points on element
        for (r = 0; r < q.n; r++) {
            uquad[r] = eval(unode,q.xi[r],q.eta[r]);
            xx = aloc[en[0]].x + dx1 * q.xi[r] + dx2 * q.eta[r];
            yy = aloc[en[0]].y + dy1 * q.xi[r] + dy2 * q.eta[r];
            aquad[r] = user->a_fcn(uquad[r],xx,yy);
        }
        // generate 3x3 element stiffness matrix (3x3 is max size but may be smaller)
        cr = 0;  cv = 0;  // cr = count rows; cv = value counter for all entries
        for (l = 0; l < 3; l++) {
            if (abf[en[l]] != 2) {
                row[cr++] = en[l];
                for (m = 0; m < 3; m++) {
                    if (abf[en[m]] != 2) {
                        sum = 0.0;
                        for (r = 0; r < q.n; r++) {
                            sum += q.w[r] * aquad[r] * InnerProd(gradpsi[l],gradpsi[m]);
                        }
                        v[cv++] = fabs(detJ) * sum;
                    }
                }
            }
        }
        // insert element stiffness matrix
        ierr = MatSetValues(P,cr,row,cr,row,v,ADD_VALUES); CHKERRQ(ierr);
    }
    ierr = ISRestoreIndices(user->mesh->e,&ae); CHKERRQ(ierr);
    ierr = ISRestoreIndices(user->mesh->bf,&abf); CHKERRQ(ierr);
    ierr = VecRestoreArrayRead(u,&au); CHKERRQ(ierr);
    ierr = UMRestoreNodeCoordArrayRead(user->mesh,&aloc); CHKERRQ(ierr);

    ierr = MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    ierr = MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    if (A != P) {
        ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
        ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    }
    ierr = MatSetOption(P,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE); CHKERRQ(ierr);
    PetscLogStagePop();  //STRIP
    return 0;
}

/* Note that nnz[n] is the number of nonzeros in row n.  It equals one for
Dirichlet rows, one more than the number of incident triangles for an
interior point, and two more than the number of incident triangles for
Neumann boundary nodes. */
//STARTPREALLOC
PetscErrorCode Preallocation(Mat J, unfemCtx *user) {
    PetscErrorCode ierr;
    const int    *ae, *abf, *en;
    int          *nnz, n, k, l;

    nnz = (int *)malloc(sizeof(int)*(user->mesh->N));
    ierr = ISGetIndices(user->mesh->bf,&abf); CHKERRQ(ierr);
    for (n = 0; n < user->mesh->N; n++)
        nnz[n] = (abf[n] == 1) ? 2 : 1;
    ierr = ISGetIndices(user->mesh->e,&ae); CHKERRQ(ierr);
    for (k = 0; k < user->mesh->K; k++) {
        en = ae + 3*k;  // en[0], en[1], en[2] are nodes of element k
        for (l = 0; l < 3; l++)
            if (abf[en[l]] != 2)
                nnz[en[l]] += 1;
    }
    ierr = ISRestoreIndices(user->mesh->e,&ae); CHKERRQ(ierr);
    ierr = ISRestoreIndices(user->mesh->bf,&abf); CHKERRQ(ierr);
    ierr = MatSeqAIJSetPreallocation(J,-1,nnz); CHKERRQ(ierr);
    free(nnz);
    return 0;
}
//ENDPREALLOC