const PluginPanelItem * operator [](size_t index) const {
			int m_ppiSize = psi().Control(m_hndl, FCTL_GETPANELITEM, index, nullptr);
			if (WinMem::Realloc(m_ppi, m_ppiSize)) {
				psi().Control(m_hndl, FCTL_GETPANELITEM, index, (LONG_PTR)m_ppi);
			}
			return m_ppi;
		}

		PCWSTR get_current_directory() const {
			int size = psi().Control(m_hndl, FCTL_GETPANELDIR, 0, nullptr);
			if (WinMem::Realloc(m_dir, size * sizeof(WCHAR), 0) && psi().Control(m_hndl, FCTL_GETPANELDIR, size, (LONG_PTR)m_dir)) {
				return	m_dir;
			}
			return L"";
		}

		const PluginPanelItem * get_selected(size_t index) const {
			int m_ppiSize = psi().Control(m_hndl, FCTL_GETSELECTEDPANELITEM, index, nullptr);
			if (WinMem::Realloc(m_ppi, m_ppiSize)) {
				psi().Control(m_hndl, FCTL_GETSELECTEDPANELITEM, index, (LONG_PTR)m_ppi);
			}
			return m_ppi;
		}

  // Determine the MA filter coefficients from the AR coefficients by
  // equating coefficients in the polynomial phi(z)*psi(z) = 1.
  // phi(z) = 1 - phi_1 z - phi_2 z^2 - ... - phi_p z^p.  This implies
  // that psi[0] = 1.  The coefficient of z^n is
  //
  // psi[n] - psi[n-1]*phi_1 - psi[n-2]*phi_2 - ... = 0
  // implying
  // psi[n] = psi[n-1]*phi_1 + ...
  //
  // The preceding math is unit-offset for phi, but zero-offset for
  // psi.
  void ArModel::set_filter_coefficients()const{
    if (filter_coefficients_current_) return;
    const Vec &phi(this->phi());
    int p = phi.size();

    filter_coefficients_.resize(2);
    filter_coefficients_[0] = 1.0;
    if(phi.empty()) return;
    filter_coefficients_[1] = phi[0];
    bool done = false;

    for(int n = 2; ; ++n){
      if (n <= phi.size()) {
        ConstVectorView phi_view(phi, 0, n);
        ConstVectorView psi(filter_coefficients_, 0, n);
        double value = phi_view.dot(psi.reverse());
        filter_coefficients_.push_back(value);
      } else {
        ConstVectorView psi(filter_coefficients_, n-p, p);
        double value = phi.dot(psi.reverse());
        filter_coefficients_.push_back(value);
        ConstVectorView psi_tail(filter_coefficients_, n-p, p);
        // You're done when the last p elements of the vector are all
        // small.
        done = psi_tail.abs_norm() < 1e-6;
      }
      if(done) break;
    }
    filter_coefficients_current_ = true;
  }

int main()
{
    double x;

    printf("  x       %-22s %-22s %-22s\n",
        "psi(x)", "psi'(x)", "psi''(x)");
    for (x = -5.5; x <= 0.5; x++)
        printf("%4.1f  % -22.15g % -22.15g % -22.15g\n",
            x, psi(x), polygamma(1, x), polygamma(2, x));
    for (x = 1; x <= 5; x++)
        printf("%4.1f  % -22.15g % -22.15g % -22.15g\n",
            x, psi(x), polygamma(1, x), polygamma(2, x));
    for (x = 10; x <= 40; x += 10)
        printf("%4.1f  % -22.15g % -22.15g % -22.15g\n",
            x, psi(x), polygamma(1, x), polygamma(2, x));
    printf("\n  x       %-22s %-22s %-22s\n",
        "psi3(x)", "psi4(x)", "psi5(x)");
    for (x = -5.5; x <= 0.5; x++)
        printf("%4.1f  % -22.15g % -22.15g % -22.15g\n",
            x, polygamma(3, x), polygamma(4, x), polygamma(5, x));
    for (x = 1; x <= 5; x++)
        printf("%4.1f  % -22.15g % -22.15g % -22.15g\n",
            x, polygamma(3, x), polygamma(4, x), polygamma(5, x));
    for (x = 10; x <= 40; x += 10)
        printf("%4.1f  % -22.15g % -22.15g % -22.15g\n",
            x, polygamma(3, x), polygamma(4, x), polygamma(5, x));
    return EXIT_SUCCESS;
}

int
Substitutions::applyTo(Term &t) const
{
    ListIterator<Substitution *> psi(*psubs);
    for ( ; !psi.done(); psi++)
    {
        if (psi()->applyTo(t) != OK)
            return(NOTOK);
    }
    return(OK);
}

		const PluginPanelItem* FacadeImpl::get_selected(size_t index) const
		{
			size_t m_ppiSize = psi().PanelControl(m_hndl, FCTL_GETSELECTEDPANELITEM, index, nullptr);
			m_ppi = static_cast<decltype(m_ppi)>(HostRealloc(heap_type, m_ppi, m_ppiSize));

			FarGetPluginPanelItem gpi = {sizeof(gpi), m_ppiSize, m_ppi};
			psi().PanelControl(m_hndl, FCTL_GETSELECTEDPANELITEM, index, &gpi);

			LogTrace(L"[%Iu] -> %p\n", index, m_ppi);
			return m_ppi;
		}

void find_most_violated_constraint(SVECTOR **fydelta, double *rhs, 
				   EXAMPLE *ex, SVECTOR *fycached, long n, 
				   STRUCTMODEL *sm, STRUCT_LEARN_PARM *sparm,
				   double *rt_viol, double *rt_psi, 
				   long *argmax_count)
     /* returns fydelta=fy-fybar and rhs scalar value that correspond
	to the most violated constraint for example ex */
{
  double      rt2=0;
  LABEL       ybar;
  SVECTOR     *fybar, *fy;
  double      factor,lossval;

  if(struct_verbosity>=2) rt2=get_runtime();
  (*argmax_count)++;
  if(sparm->loss_type == SLACK_RESCALING) 
    ybar=find_most_violated_constraint_slackrescaling(ex->x,ex->y,sm,sparm);
  else
    ybar=find_most_violated_constraint_marginrescaling(ex->x,ex->y,sm,sparm);
  if(struct_verbosity>=2) (*rt_viol)+=MAX(get_runtime()-rt2,0);
  
  if(empty_label(ybar)) {
    printf("ERROR: empty label was returned for example\n");
    /* exit(1); */
    /* continue; */
  }
  
  /**** get psi(x,y) and psi(x,ybar) ****/
  if(struct_verbosity>=2) rt2=get_runtime();
  if(fycached)
    fy=copy_svector(fycached); 
  else 
    fy=psi(ex->x,ex->y,sm,sparm);
  fybar=psi(ex->x,ybar,sm,sparm);
  if(struct_verbosity>=2) (*rt_psi)+=MAX(get_runtime()-rt2,0);
  lossval=loss(ex->y,ybar,sparm);
  free_label(ybar);
  
  /**** scale feature vector and margin by loss ****/
  if(sparm->loss_type == SLACK_RESCALING)
    factor=lossval/n;
  else                 /* do not rescale vector for */
    factor=1.0/n;      /* margin rescaling loss type */
  mult_svector_list(fy,factor);
  mult_svector_list(fybar,-factor);
  append_svector_list(fybar,fy);   /* compute fy-fybar */

  (*fydelta)=fybar;
  (*rhs)=lossval/n;
}

LABEL       find_most_violated_constraint_slackrescaling(PATTERNX x, LABEL y, 
						     STRUCTMODEL *sm, 
						     STRUCT_LEARN_PARM *sparm)
{
  /* Finds the label ybar for pattern x that that is responsible for
     the most violated constraint for the slack rescaling
     formulation. It has to take into account the scoring function in
     sm, especially the weights sm.w, as well as the loss
     function. The weights in sm.w correspond to the features defined
     by psi() and range from index 1 to index sm->sizePsi. Most simple
     is the case of the zero/one loss function. For the zero/one loss,
     this function should return the highest scoring label ybar, if
     ybar is unequal y; if it is equal to the correct label y, then
     the function shall return the second highest scoring label. If
     the function cannot find a label, it shall return an empty label
     as recognized by the function empty_label(y). */
  LABEL ybar;
  DOC doc;
  long classlabel, bestclass=-1, first=1;
  double score, score_y, score_ybar, bestscore=-1;

  /* NOTE: This function could be made much more efficient by not
     always computing a new PSI vector. */
  doc = *(x.doc);
  doc.fvec = psi(x,y,sm,sparm);
  score_y = classify_example(sm->svm_model,&doc);
  free_svector(doc.fvec);

  ybar.scores = NULL;
  ybar.num_classes = sparm->num_classes;
  for(classlabel=1; classlabel<=sparm->num_classes; classlabel++) {
    ybar.classlabel = classlabel;
    doc.fvec=psi(x,ybar,sm,sparm);
    score_ybar=classify_example(sm->svm_model,&doc);
    free_svector(doc.fvec);
    score=loss(y,ybar,sparm,x.doc->fvec)*(1.0-score_y+score_ybar);
    if((bestscore<score)  || (first)) {
      bestscore=score;
      bestclass = classlabel;
      first=0;
    }
  }
  if(bestclass == -1) 
    printf("ERROR: Only one class\n");
  ybar.classlabel = bestclass;
  if(struct_verbosity>=3)
    printf("[%ld:%.2f] ",bestclass,bestscore);
  return(ybar);
}

vector horizontalVelocityField::streamfunctionAt
(
        const point& p,
        const Time& t
) const
{
    const vector unitNormal(0, -1, 0);
    const dimensionedScalar z("z", dimLength, p.z());

    dimensionedScalar psi("psi", cmptMultiply(dimVelocity, dimLength), scalar(0));
    if (z.value() <= z1.value())
    {
        // psi is zero
    }
    else if (z.value() <= z2.value())
    {
        psi = -0.5*u0*(z - z1 - (z2-z1)/M_PI*Foam::sin(M_PI*(z-z1)/(z2-z1)));
    }
    else 
    {
        psi = -0.5*u0*(2*z - z2 - z1);
    }

    return unitNormal * psi.value();
}

	inline void ebox_code(DWORD err, PCWSTR line) {
		ustring	title(L"Error: ");
		title += as_str(err);
		::SetLastError(err);
		PCWSTR Msg[] = {title.c_str(), line, L"OK", };
		psi().Message(get_plugin_guid(), FMSG_WARNING | FMSG_ERRORTYPE, nullptr, Msg, sizeofa(Msg), 1);
	}

	bool Settings_t::set(PCWSTR name, PCWSTR value, FARSETTINGS_SUBFOLDERS root)
	{
		LogTrace();
		FarSettingsItem item = {sizeof(item), root, name, FST_STRING};
		item.String = value;
		return psi().SettingsControl(m_hndl, SCTL_SET, 0, &item);
	}

		Panel(const HANDLE aPlugin, int cmd = FCTL_GETPANELINFO):
			m_hndl(aPlugin),
			m_ppi(nullptr),
			m_dir(nullptr) {
			WinMem::Zero(m_pi);
			m_Result = psi().Control(m_hndl, cmd, 0, (LONG_PTR)&m_pi);
		}

void test_staggered() {
  mdp << "START TESTING STAGGERED ACTIONS\n";

  int box[]={64,6,6,6}, nc=3;
  generic_lattice lattice(4,box,default_partitioning<0>,
			  torus_topology, 0, 3);
  gauge_field U(lattice,nc);
  gauge_field V(lattice,nc);
  staggered_field psi(lattice, nc);
  staggered_field chi1(lattice, nc);
  staggered_field chi2(lattice, nc);
  coefficients coeff;
  coeff["mass"]=1.0;
  double t0, t1;
  inversion_stats stats;
  set_hot(U);
  set_random(psi);

  mdp << "ATTENTION: need to adjust asqtad coefficnets\n";

  default_staggered_action=StaggeredAsqtadActionFast::mul_Q;

  default_staggered_inverter=MinimumResidueInverter<staggered_field,gauge_field>;
  t0=mpi.time();
  stats=mul_invQ(chi2,psi,U,coeff);
  t1=(mpi.time()-t0)/lattice.nvol_gl/stats.steps;
  cout << "Staggered Min Res TIME=" << t1 << endl;

  default_staggered_inverter=BiConjugateGradientStabilizedInverter<staggered_field,gauge_field>;
  t0=mpi.time();
  stats=mul_invQ(chi2,psi,U,coeff);
  t1=(mpi.time()-t0)/lattice.nvol_gl/stats.steps;
  cout << "Staggered BiCGStab TIME=" << t1 << endl;

  default_staggered_inverter=StaggeredBiCGUML::inverter;
  t0=mpi.time();
  stats=mul_invQ(chi2,psi,U,coeff);
  t1=(mpi.time()-t0)/lattice.nvol_gl/stats.steps;
  cout << "Staggered SSE BiCGStabUML TIME=" << t1 << endl;

  default_staggered_action=StaggeredAsqtadActionSSE2::mul_Q;

  default_staggered_inverter=MinimumResidueInverter<staggered_field,gauge_field>;
  t0=mpi.time();
  stats=mul_invQ(chi2,psi,U,coeff);
  t1=(mpi.time()-t0)/lattice.nvol_gl/stats.steps;
  cout << "Staggered SSE Min Res TIME=" << t1 << endl;

  default_staggered_inverter=BiConjugateGradientStabilizedInverter<staggered_field,gauge_field>;
  t0=mpi.time();
  stats=mul_invQ(chi2,psi,U,coeff);
  t1=(mpi.time()-t0)/lattice.nvol_gl/stats.steps;
  cout << "Staggered SSE BiCGStab TIME=" << t1 << endl;

  default_staggered_inverter=StaggeredBiCGUML::inverter;
  t0=mpi.time();
  stats=mul_invQ(chi2,psi,U,coeff);
  t1=(mpi.time()-t0)/lattice.nvol_gl/stats.steps;
  cout << "Staggered SSE BiCGStabUML TIME=" << t1 << endl;
}

double
xxliked (int m, double a, double logsum, double lsum)
// first deriv wrt a
{
  double p, yl = 0.0;
  double ym, x, var, vard;
  int j;

  p = 0.5 * (double) (m + 1);
  ym = (double) m;
  var = lsum / (2.0 * a * ym);
  vard = -var / a;

  yl = -ym * log (2.0);
  for (j = 1; j <= m; ++j)
    {
      x = a - 0.5 * (double) (m - j);
      if (x < 0.0)
        return 100.0;
      yl -= psi (x);
    }
// so far this is log (C_L)    normalizing constant
  yl -= ym * log (var);
  yl -= (ym * a / var) * vard;
  yl += logsum;
  yl -= ym;  // plugging in var 

  return yl;

}

bool DGKTSD::evaluateA(const DKTS& a, const DKTS& xi, 
                       const DKTS& v, DKTS& w,
                       const LongReal& alpha, const LongReal& beta) const
 {
    bool value = true;
    LongInt d = DGKTSD::d();
    LongInt k = DGKTSD::k();

    w().setNull();

    for(LongInt mu1=0; mu1<d; mu1++)
     {
        for(LongInt j1=0; j1<k; j1++)
         {
            DKTVector& w_jmu = w(j1, mu1);
            const LongReal& preFactor = W(j1, mu1);

            // The Matrix A
            for(LongInt j2=0; j2<k; j2++)
             {
                w_jmu.update(psi(j1, j2, mu1, mu1)*W(j2, mu1), v(j2, mu1));
             }

            w_jmu *= preFactor*beta;
         }// End for(LongInt j1=0; j1<k; j1++)
     }// End for(LongInt mu1=0; mu1<d; mu1++)    

    w().update(alpha, v());

   return value;
 }

 double SBShapelet::SBShapeletImpl::xValue(const Position<double>& p) const
 {
     LVector psi(_bvec.getOrder());
     psi.fillBasis(p.x/_sigma, p.y/_sigma, _sigma);
     double xval = _bvec.dot(psi);
     return xval;
 }

void test_psi() {
    int i;
    for(i=1; i<1000; i++) {
        printf("%f, ",psi((double)i));
        }

    }

		intptr_t insert_string(ssize_t y, PCWSTR str, int size, PCWSTR eol)
		{
			if (set_position(y) && psi().EditorControl(-1, ECTL_INSERTSTRING, 0, 0)) {
				return set_string(y, str, size, eol);
			}
			return 0;
		}

		ssize_t get_string(ssize_t y, PCWSTR & str)
		{
			EditorGetString egs = {sizeof(egs), y};
			psi().EditorControl(-1, ECTL_GETSTRING, 0, &egs);
			str = egs.StringText;
			return egs.StringLength;
		}