/* ----------------------------------------------------------------------------

Copyright (C) 2009.

A. Ronald Gallant
Post Office Box 659
Chapel Hill NC 27514-0659
USA

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.

This program 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 General Public License for more details.

You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

-----------------------------------------------------------------------------*/

#include "libsmm.h"
#include "emm.h"

using namespace std;
using namespace scl;
using namespace emm;
using namespace libsmm;

namespace emm {
  
  vector<string> 
  habit_usrmod_variables::get_habit_usrmod_variables(realmat& mv)
  {
    vector<string> names(17);
    names[0] = "log_consumption_growth";
    mv = log_consumption_growth;
    names[1] = "log_dividend_growth";
    mv = cbind(mv,log_dividend_growth);
    names[2] = "log_surplus_ratio";
    mv = cbind(mv,log_surplus_ratio);
    names[3] = "consumption_growth";
    mv = cbind(mv,consumption_growth);
    names[4] = "dividend_growth";
    mv = cbind(mv,dividend_growth);
    names[5] = "consumption";
    mv = cbind(mv,consumption);
    names[6] = "dividend";
    mv = cbind(mv,dividend);
    names[7] = "surplus_ratio";
    mv = cbind(mv,surplus_ratio);
    names[8] = "sensitivity_function";
    mv = cbind(mv,sensitivity_function);
    names[9] = "marginal_rate_of_substitution";
    mv = cbind(mv,marginal_rate_of_substitution);
    names[10] = "price_dividend_ratio";
    mv = cbind(mv,price_dividend_ratio);
    names[11] = "bond_price"; 
    mv = cbind(mv,bond_price);
    names[12] = "gross_stock_return";
    mv=cbind(mv,gross_stock_return);
    names[13] = "gross_risk_free_rate";
    mv=cbind(mv,gross_risk_free_rate);
    names[14] = "geometric_stock_return";
    mv=cbind(mv,geometric_stock_return);
    names[15] = "geometric_risk_free_rate";
    mv=cbind(mv,geometric_risk_free_rate);
    realmat log_price_dividend_ratio = price_dividend_ratio;
    for (INTEGER i=1; i<=3; ++i) log_price_dividend_ratio[i] = 0.0;
    for (INTEGER i=4; i<=nt; ++i) 
      log_price_dividend_ratio[i] = log(log_price_dividend_ratio[i]);
    names[16] = "log_price_dividend_ratio";
    mv = cbind(mv,log_price_dividend_ratio);
    return names;
  }
  
  habit_usrmod::habit_usrmod(const scl::realmat& dat, INTEGER len_mod_parm,
      INTEGER len_mod_func, const std::vector<std::string>& pfvec,
      const std::vector<std::string>& alvec, std::ostream& detail)
  : mp(), mv(120,60000,120), data(dat), 
    bootstrap_seed(770116), simulation_seed(740726)
  {
    if (alvec.size() < 3) error("Error, habit_usrmod constructor, bad alvec");

    vector<string>::const_iterator alvec_ptr = alvec.begin();

    INTEGER n0 = atoi((++alvec_ptr)->substr(0,12).c_str());
    INTEGER n  = atoi((++alvec_ptr)->substr(0,12).c_str());
    INTEGER n1 = atoi((++alvec_ptr)->substr(0,12).c_str());

    mv = habit_usrmod_variables(n0,n,n1);
 
    realmat sgrid(22,1), pdval(22,1);  // CC policy function
  
    sgrid[1] = -5.20514;  pdval[1] = 4.35107;
    sgrid[2] = -4.96183;  pdval[2] = 4.44397;
    sgrid[3] = -4.76405;  pdval[3] = 4.52012;
    sgrid[4] = -4.74214;  pdval[4] = 4.5287;
    sgrid[5] = -4.52244;  pdval[5] = 4.61509;
    sgrid[6] = -4.27914;  pdval[6] = 4.7138;
    sgrid[7] = -4.07091;  pdval[7] = 4.80184;
    sgrid[8] = -3.66544;  pdval[8] = 4.9827;
    sgrid[9] = -3.37776;  pdval[9] = 5.11901;
    sgrid[10] = -3.15462;  pdval[10] = 5.23134;
    sgrid[11] = -2.97229;  pdval[11] = 5.32856;
    sgrid[12] = -2.86453;  pdval[12] = 5.38912;
    sgrid[13] = -2.81814;  pdval[13] = 5.41611;
    sgrid[14] = -2.68461;  pdval[14] = 5.49662;
    sgrid[15] = -2.56683;  pdval[15] = 5.57113;
    sgrid[16] = -2.46147;  pdval[16] = 5.64058;
    sgrid[17] = -2.40616;  pdval[17] = 5.67804;
    sgrid[18] = -2.39616;  pdval[18] = 5.6849;
    sgrid[19] = -2.38616;  pdval[19] = 5.69179;
    sgrid[20] = -2.37616;  pdval[20] = 5.69873;
    sgrid[21] = -2.37191;  pdval[21] = 5.7017;
    sgrid[22] = -2.36616;  pdval[22] = 5.70575;
  
    CCpdval.update(sgrid,pdval);

    if (n_parms != len_mod_parm) {
      error("Error, usrmod, constructor, len_mod_parm is wrong in parmfile");
    }

    if (n_stats != len_mod_func) {
      error("Error, usrmod, constructor, len_mod_parm is wrong in parmfile");
    }

    if (n_datum != dat.nrow()) {
      error("Error, usrmod, constructor, M is wrong in parmfile");
    }

  }
  
  void habit_usrmod::make_state(INT_32BIT& seed)
  {
    const INT_32BIT start_seed = 100542;
  
    REAL del_c;
    REAL del_d;
    
    REAL cg;
    REAL dg;
    
    REAL s; 
    REAL S_bar = mp.sig*sqrt(mp.gamma/(1.0-mp.phi));
    REAL s_bar = log(S_bar);
    REAL s_max = s_bar + (0.5)*(1.0 - pow(S_bar,2));
    REAL s_lag = s_bar;
    REAL del_s;
  
    REAL lam; 
    REAL lam_lag = 0.0;
  
    REAL mrs;
    
    realmat e(2,1);
    
    INT_32BIT state_seed = start_seed;
  
    REAL C = REAL_MIN;
    REAL D = REAL_MIN;

    for (INTEGER t=2; t <=mv.nt; t++) {
  
      if (t == mv.n0+1) state_seed = seed;
  
      e[1] = unsk(state_seed);
      e[2] = unsk(state_seed);
  
      e = mp.R*e;
   
      del_c = mp.g + e[1];
      del_d = mp.g + e[2];
  
      cg = exp(del_c);
      dg = exp(del_d);
    
      C *= cg;
      D *= dg;
  
      mv.log_consumption_growth[t] = del_c;
      mv.consumption_growth[t] = cg;
      mv.consumption[t] = C;
      
      mv.log_dividend_growth[t] = del_d;
      mv.dividend_growth[t] = dg;
      mv.dividend[t] = D;
      
      s = (1.0 - mp.phi)*s_bar + mp.phi*s_lag + lam_lag*e[1];
      del_s = s - s_lag;
      s_lag = s;
      
      mv.log_surplus_ratio[t] = s;
      mv.surplus_ratio[t] = exp(s);
  
      lam = (s < s_max) ? sqrt(1.0 - 2.0*(s-s_bar))/S_bar - 1.0 : 0.0;
      lam_lag = lam;
      
      mv.sensitivity_function[t] = lam;
  
      mrs = mp.delta*exp(-mp.gamma*(del_s + del_c));
  
      mv.marginal_rate_of_substitution[t] = mrs;
    }
    seed = state_seed;
  }
  
  bool habit_usrmod::make_sim(INT_32BIT& seed, realmat& sim)
  {
    realmat theta(mp.p,1);
    this -> get_rho(theta);
      
    this->make_state(seed);
  
    REAL smax = -REAL_MAX;
    for (INTEGER t=mv.n0+1; t<=mv.n0+mv.n; ++t) {
      smax = smax < mv.log_surplus_ratio[t] ? mv.log_surplus_ratio[t] : smax;
    }
  
    REAL S_bar = mp.sig*sqrt(mp.gamma/(1.0-mp.phi));
    REAL s_bar = log(S_bar);
    REAL s_max = s_bar + (0.5)*(1.0 - pow(S_bar,2));
    REAL S_max = exp(s_max);
  
    smax = smax < s_max ? s_max : smax;
  
    REAL top = exp(smax);
  
    const INTEGER intervals = 10;        // the following construction of
                                         // sgrid replicates that in CC
    REAL increment = top/(intervals+1);
    
    bool no_S_bar = true;
    bool no_S_max = true;;
    realmat sgrid;
    for (INTEGER i=1; i<=intervals; ++i) {  // equal increments in S
      REAL Si = i*increment;
      if (Si == S_bar) no_S_bar = false;
      if (Si == S_max) no_S_max = false;
      sgrid.push_back(log(Si));
    }
    if (no_S_bar) sgrid.push_back(s_bar);  // make sure these two in grid
    if (no_S_max) sgrid.push_back(s_max);
    sgrid.push_back(smax - 0.01);          // extra points added near smax
    sgrid.push_back(smax - 0.02);
    sgrid.push_back(smax - 0.03);
    sgrid.push_back(smax - 0.04);
  
    sgrid.sort();
    REAL smin = sgrid[1];
  
    const INTEGER npoints = 5;    // adding quadrature points to sgrid
                                  // is a departure from CC
    realmat x, w;
    if( hquad(npoints,x,w) ) error("Error, usrmod, gen_sim, hquad failed");
    
    const REAL root2 = sqrt(2.0);              // 1.4142135623730951
    const REAL rootpi = sqrt(4.0*atan(1.0));   // 1.7724538509055161
    
    realmat abscissae(npoints,1);
    realmat weights(npoints,1);
    for (INTEGER i=1; i<=npoints; ++i) {
      abscissae[i] = root2*mp.sig*x[i];
      weights[i] = w[i]/rootpi;
    }
  
    for (INTEGER j=1; j<=npoints; ++j) {
      REAL s0 = smin;
      REAL lam = (s0 < s_max) ? sqrt(1.0 - 2.0*(s0-s_bar))/S_bar - 1.0 : 0.0;
      REAL e1 = abscissae[j];
      REAL s1 = (1.0 - mp.phi)*s_bar + mp.phi*s0 + lam*e1;
      sgrid.push_back(s1);
      s0 = smax;
      lam = (s0 < s_max) ? sqrt(1.0 - 2.0*(s0-s_bar))/S_bar - 1.0 : 0.0;
      e1 = abscissae[j];
      s1 = (1.0 - mp.phi)*s_bar + mp.phi*s0 + lam*e1;
      sgrid.push_back(s1);
    }
    
    sgrid.sort();
  
    INTEGER ngrid = sgrid.size();          // delete points that are too close
    intvec idx = seq(1,ngrid);
    for (INTEGER i=1; i<ngrid; ++i) {
      if (sgrid[i+1] - sgrid[i] < 0.001) idx[i+1] = -1;
    }
  
    sgrid = sgrid(idx,"");
    ngrid = sgrid.size();
  
    intvec tstidx = seq(1,ngrid);       // only check convergence between
    INTEGER mid = (ngrid + 1)/2;        // smin and smax 
    for (INTEGER i=1; i<=ngrid; ++i) {   
      if (sgrid[i] <= smin || smax <= sgrid[i]) tstidx[i] = mid; 
    }
  
    realmat pdval(ngrid,1);
  
    for (INTEGER i=1; i<=ngrid; ++i) {
      pdval[i] = CCpdval(sgrid[i]);
    }
    
    linear_interpolater v(sgrid,pdval);  
    
    // v maps log_surplus_ratio to log_price_dividend_ratio, i.e.,
    // log_price_dividend_ratio = v(log_surplus_ratio)
    
    const INTEGER maxiter = 5000;
    const REAL tol = 1.0e-4;
    bool converge = false;
  
    realmat PD0(ngrid,1);
    realmat PD0_lag(ngrid,1);
  
    for (INTEGER i=1; i<=ngrid; ++i) {
      PD0_lag[i] = exp(v(sgrid[i]));
    }
  
    REAL Ewv_fixed_part = exp(0.5*(1.0 - pow(mp.rho,2))*pow(mp.sig_w,2));
    REAL G = exp(mp.g);
  
    for (INTEGER iter=1; iter<=maxiter; ++iter) {
      for (INTEGER i=1; i<=ngrid; ++i) {
        REAL s0 = sgrid[i];
        REAL lam = (s0 < s_max) ? sqrt(1.0 - 2.0*(s0-s_bar))/S_bar - 1.0 : 0.0;
        PD0[i] = 0.0;
        for (INTEGER j=1; j<=npoints; ++j) {
          REAL e1 = abscissae[j];
          REAL s1 = (1.0 - mp.phi)*s_bar + mp.phi*s0 + lam*e1;
          REAL del_s = s1 - s0;
          REAL del_c = mp.g + e1;
          REAL M1 = mp.delta*exp(-mp.gamma*(del_s + del_c));
          REAL PD1 = exp(v(s1));
          REAL Ewv = Ewv_fixed_part*exp(mp.rho*(mp.sig_w/mp.sig)*e1);
          PD0[i] += G*M1*Ewv*(1.0 + PD1)*weights[j];
        }
      }
  
      for (INTEGER i=1; i<=ngrid; ++i) {
        pdval[i] = log(PD0[i]);
      }
      v.update(sgrid,pdval);
  
      converge = true;
      for (INTEGER i=1; i<=ngrid; ++i) {
        INTEGER j = tstidx[i];
        if (fabs(PD0[j] - PD0_lag[j]) > fabs(PD0_lag[j] + tol)*tol) {
          converge = false;
          break;
        }
      }
  
      for (INTEGER i=1; i<=ngrid; ++i) {
        PD0_lag[i] = PD0[i];
      }
      
      if (converge) break;
    }
  
    if (!converge) {
      warn("Warning, convergence failed at this parameter value");
      cerr << "\t g     " << theta[1] << '\n';
      cerr << "\t R[1]  " << theta[2] << '\n';
      cerr << "\t R[2]  " << theta[3] << '\n';
      cerr << "\t R[3]  " << theta[4] << '\n';
      cerr << "\t phi   " << theta[5] << '\n';
      cerr << "\t delta " << theta[6] << '\n';
      cerr << "\t gamma " << theta[7] << '\n';
      cerr << "\t mudc  " << theta[8] << '\n';
    }
  
    realmat bval(ngrid,1,-0.00119638);    // fill correspondes to 
                                          // CC parameter values
    linear_interpolater b(sgrid,bval);
    
    // b maps log_surplus_ratio to the log_bond_price, i.e.,
    // log_bond_price = B(log_surplus_ratio)
  
    realmat B0(ngrid,1);
    
    for (INTEGER i=1; i<=ngrid; ++i) {
      REAL s0 = sgrid[i];
      REAL lam = (s0 < s_max) ? sqrt(1.0 - 2.0*(s0-s_bar))/S_bar - 1.0 : 0.0;
      B0[i] = 0.0;
      for (INTEGER j=1; j<=npoints; ++j) {
        REAL e1 = abscissae[j];
        REAL s1 = (1.0 - mp.phi)*s_bar + mp.phi*s0 + lam*e1;
        REAL del_s = s1 - s0;
        REAL del_c = mp.g + e1;
        REAL M1 = mp.delta*exp(-mp.gamma*(del_s + del_c));
        B0[i] += M1*weights[j];
      }
    }
  
    for (INTEGER i=1; i<=ngrid; ++i) {
      bval[i] = log(B0[i]);
    }
    
    b.update(sgrid,bval);
  
    for (INTEGER t=2; t <= mv.nt; ++t) {
      mv.price_dividend_ratio[t] = exp(v(mv.log_surplus_ratio[t]));
      mv.bond_price[t] = exp(b(mv.log_surplus_ratio[t]));
    }
  
    for (INTEGER t=3; t <= mv.nt; ++t) {
      REAL ret = 1.0 + mv.price_dividend_ratio[t];
      ret /= mv.price_dividend_ratio[t-1];
      ret *= mv.dividend_growth[t];
      mv.gross_stock_return[t] = ret;
      mv.gross_risk_free_rate[t] = (1.0)/mv.bond_price[t];
      mv.geometric_stock_return[t] = log(ret);
      mv.geometric_risk_free_rate[t] = -log(mv.bond_price[t]);
    }
  
    sim.resize(n_datum,mv.n/12);
  
    REAL r0,C0,C12,D0,D12,PD;     
    INTEGER tt = 0;
    for (INTEGER i=1; i <= mv.n/12; ++i) {
      r0 = C0 = C12 = D0 = D12 = 0.0;
      for (INTEGER j=1; j <= 12; j++) {
        ++tt;
        r0  += mv.geometric_stock_return[mv.n0 + tt];
        C0  += mv.consumption[mv.n0 + tt];
        C12 += mv.consumption[mv.n0 + tt - 12];
        D0  += mv.dividend[mv.n0 + tt];
        D12 += mv.dividend[mv.n0 + tt - 12];
      }
      PD = mv.price_dividend_ratio[mv.n0 + tt]*mv.dividend[mv.n0 + tt]/D0;
      sim(1,i) = log(D0) - log(C0) + mp.mudc;
      sim(2,i) = log(C0) - log(C12);
      sim(3,i) = log(PD);
      sim(4,i) = r0;
    }
  
    return converge; 
  }
  
  bool habit_usrmod::gen_sim(realmat& sim)
  {
    INT_32BIT seed = simulation_seed;
    return make_sim(seed, sim);
  }
  
  bool habit_usrmod::gen_bootstrap(vector<realmat>& bs)
  {
    INTEGER lrho = mp.p;
    vector<realmat>::size_type len_vec = 2*lrho+1;
  
    if (len_vec != bs.size()) {
      bs.resize(len_vec);
    }
  
    INTEGER n_dat = data.ncol();
    INTEGER r_dat = data.nrow();
    INTEGER len = len_vec*n_dat;
  
    realmat sim;
    bool converge =  make_sim(bootstrap_seed, sim);   
    if (!converge) return false;
  
    if (len < n_dat) {
      string msg = "Error, gen_bootstrap, simulation length must be at least";
      msg += fmt('d',7,len).get_ostr();
      error(msg);
    }
  
    realmat dat(r_dat,n_dat);
    INTEGER s = 0;
    for (vector<realmat>::size_type k=0; k<len_vec; ++k) {
      for (INTEGER j=1; j<=n_dat; ++j) {
        for (INTEGER i=1; i<=r_dat; ++i) {
          dat(i,j) = sim(i,s+j);
        }
      }
      bs[k] = dat;
      s += n_dat;
    }
  
    return true;
  }
  
  bool habit_usrmod::gen_sim(realmat& sim, realmat& stats)
  {
    bool success = this->gen_sim(sim);
  
    if (success) {
    
      REAL mean_rf = 0.0;
      REAL mean_sr = 0.0;
  
      for (INTEGER t=mv.n0+1; t <= mv.n0+mv.n; t++) {
        mean_rf += mv.geometric_risk_free_rate[t];
        mean_sr += mv.geometric_stock_return[t];
      }
  
      mean_rf /= REAL(mv.n);
      mean_sr /= REAL(mv.n);
    
      REAL var_rf = 0.0;
      REAL var_sr = 0.0;
      
      for (INTEGER t=mv.n0+1; t <= mv.n0+mv.n; t++) {
        var_rf += pow(mv.geometric_risk_free_rate[t]-mean_rf, 2);
        var_sr += pow(mv.geometric_stock_return[t]-mean_sr, 2);
      }
  
      var_rf /= REAL(mv.n);
      var_sr /= REAL(mv.n);
    
      stats.resize(n_stats,1);
  
      stats.check(1) = mean_rf;
      stats.check(2) = var_rf;
      stats.check(3) = mean_sr;
      stats.check(4) = var_sr;
  
    }
    return success;
  }
  
  bool habit_usrmod::support(const realmat& theta) 
  {
    realmat theta_lo(n_parms,1);
    
    theta_lo[1] = -0.005;  // g 
    theta_lo[2] =  0.0;    // R(1,1)
    theta_lo[3] = -1.0;    // R(1,2)
    theta_lo[4] =  0.0;    // R(2,2)
    theta_lo[5] = -1.0;    // phi  
    theta_lo[6] =  0.9;    // delta
    theta_lo[7] =  0.25;   // gamma
    theta_lo[8] = -10.0;   // mudc 
  
    realmat theta_hi(n_parms,1);
  
    theta_hi[1] =  0.005;  // g 
    theta_hi[2] =  1.0;    // R(1,1)
    theta_hi[3] =  1.0;    // R(1,2)
    theta_hi[4] =  1.0;    // R(2,2)
    theta_hi[5] =  1.0;    // phi  
    theta_hi[6] =  1.05;   // delta
    theta_hi[7] = 10.0;    // gamma
    theta_hi[8] = 10.0;    // mudc 
  
    return ( (theta_lo < theta) && (theta < theta_hi) );
  }
  
  den_val habit_usrmod::prior(const realmat& theta,const realmat& stats)
  {
    habit_usrmod_parameters parm;
    parm.set_theta(theta);

    realmat rho(parm.p,1);
    rho[1] = parm.g;
    rho[2] = parm.sig;
    rho[3] = parm.sig_w;
    rho[4] = parm.rho;
    rho[5] = parm.phi;
    rho[6] = parm.delta;
    rho[7] = parm.gamma;
    rho[8] = parm.mudc;

    const REAL minus_log_root_two_pi = -9.1893853320467278e-01;
  
    const REAL bond = 0.000743618;             // 0.89% annualized
    const REAL sbond = (1.0/1200.0)/1.96;      // P(within 1% annualized)=0.95
  
    realmat mean(parm.p,1,0.0);
    realmat sdev(parm.p,1,100.0);

    mean[ 1] =   0.001575;      // g
    mean[ 2] =   0.00433013;    // sig
    mean[ 3] =   0.2;           // rho 
    mean[ 4] =   0.0323316;     // sig_w
    mean[ 5] =   0.988462;      // phi
    mean[ 6] =   0.990336;      // delta
    mean[ 7] =   2.0;           // gamma
    mean[ 8] =  -3.3965;        // mudc

    for (INTEGER i=1; i<=parm.p; ++i) {
      sdev[i] = fabs(mean[i])*0.1/1.96;
    }

    sdev[5] = fabs(mean[5])*0.001/1.96;    // Don't forget to report
    sdev[6] = fabs(mean[6])*0.001/1.96;    // this in the paper.

    den_val sum(true,0.0);
    
    REAL zbond = (stats[1]-bond)/sbond;
    REAL ebond = minus_log_root_two_pi - log(sbond) - 0.5*pow(zbond,2);
  
    sum += den_val(true,ebond);
  
    for (INTEGER i=1; i<=parm.p; ++i) {
      REAL z = (rho[i] - mean[i])/sdev[i];
      REAL e = minus_log_root_two_pi - log(sdev[i]) - 0.5*pow(z,2);
      sum += den_val(true,e);
    }

    return sum;
  }

  void  habit_usrmod::write_usrvar(const char* filename) 
  { 
    ofstream fout(filename);
    fout << starbox("/mean of model variables//") << '\n';
    realmat vars;
    vector<string> names = mv.get_habit_usrmod_variables(vars);
    INTEGER c = vars.ncol();
    INTEGER r = vars.nrow();
    realmat ones(1,r,1.0);
    realmat mean = ones*vars/r;
    INTEGER len = 0;
    for (INTEGER i=1; i<=c; ++i) {
      INTEGER sz = names[i-1].size();
      len = len > sz ? len :  sz;
    }
    for (INTEGER i=1; i<=c; ++i) {
      string pad(len- names[i-1].size()+2,' ');
      fout << names[i-1] << pad << mean[i] << '\n';
    }
    fout << '\n';

    fout << starbox("/mean of data//");
    r = data.nrow();
    c = data.ncol();
    ones.resize(c,1,1.0);
    mean = data*ones/c;
    fout << mean;
    fout << '\n';

    fout << starbox("/mean of annual simulation//");
    realmat sim, func;
    bool success = gen_sim(sim,func);
    if (!success) error("Error, write_usrvar, gen_sim failed");
    r = sim.nrow();
    c = sim.ncol();
    ones.resize(c,1,1.0);
    mean = sim*ones/c;
    fout << mean;
    fout << '\n';

    fout << starbox("/annual functionals//");
    r = func.nrow();
    c = func.ncol();
    ones.resize(c,1,1.0);
    mean = func*ones/c;
    fout << mean;
    return; 
  }
}