/* Absolute value versions of the above */
static int32_t MaxAbsAccumulate(CSOUND *csound, MINMAXACCUM *p)
{
IGN(csound);
uint32_t offset = p->h.insdshead->ksmps_offset;
uint32_t early = p->h.insdshead->ksmps_no_end;
uint32_t n, nsmps = CS_KSMPS;
MYFLT *out = p->accum;
MYFLT *in = p->ain;
MYFLT inabs;
if (UNLIKELY(offset)) memset(out, '\0', offset*sizeof(MYFLT));
if (UNLIKELY(early)) {
nsmps -= early;
memset(&out[nsmps], '\0', early*sizeof(MYFLT));
}
for (n=offset; n<nsmps; n++) {
inabs = FABS(in[n]);
if (UNLIKELY(inabs > out[n]))
out[n] = inabs;
}
return OK;
}
开发者ID:csound，项目名称:csound，代码行数:24，代码来源:minmax.c
示例2: CalcUnpred
// input : current spectrum in the form of power *spec and phase *phase,
// the last two earlier spectrums are at position
// 512 and 1024 of the corresponding Input-Arrays.
// Array *vocal, which can mark an FFT_Linie as harmonic
// output: current amplitude *amp and unpredictability *cw
static void
CalcUnpred (PsyModel* m,
const int MaxLine,
const float* spec,
const float* phase,
const int* vocal,
float* amp0,
float* phs0,
float* cw )
{
int n;
float amp;
float tmp;
#define amp1 ((amp0) + 512) // amp[ 512...1023] contains data of frame-1
#define amp2 ((amp0) + 1024) // amp[1024...1535] contains data of frame-2
#define phs1 ((phs0) + 512) // phs[ 512...1023] contains data of frame-1
#define phs2 ((phs0) + 1024) // phs[1024...1535] contains data of frame-2
for ( n = 0; n < MaxLine; n++ ) {
tmp = COSF ((phs0[n] = phase[n]) - 2*phs1[n] + phs2[n]); // copy phase to output-array, predict phase and calculate predictive error
amp0[n] = SQRTF (spec[n]); // calculate and set amplitude
amp = 2*amp1[n] - amp2[n]; // predict amplitude
// calculate unpredictability
cw[n] = SQRTF (spec[n] + amp * (amp - 2*amp0[n] * tmp)) / (amp0[n] + FABS(amp));
}
// postprocessing of harmonic FFT-lines (*cw is set to CVD_UNPRED)
if ( m->CVD_used && vocal != NULL ) {
for ( n = 0; n < MAX_CVD_LINE; n++, cw++, vocal++ )
if ( *vocal != 0 && *cw > CVD_UNPRED * 0.01 * *vocal )
*cw = CVD_UNPRED * 0.01 * *vocal;
}
return;
}
/* central finite difference approximation to the jacobian of func */
void FDIF_CENT_JAC_APPROX(
void (*func)(LM_REAL *p, LM_REAL *hx, int m, int n, void *adata),
/* function to differentiate */
LM_REAL *p, /* I: current parameter estimate, mx1 */
LM_REAL *hxm, /* W/O: work array for evaluating func(p-delta), nx1 */
LM_REAL *hxp, /* W/O: work array for evaluating func(p+delta), nx1 */
LM_REAL delta, /* increment for computing the jacobian */
LM_REAL *jac, /* O: array for storing approximated jacobian, nxm */
int m,
int n,
void *adata)
{
register int i, j;
LM_REAL tmp;
register LM_REAL d;
for(j=0; j<m; ++j){
/* determine d=max(1E-04*|p[j]|, delta), see HZ */
d=CNST(1E-04)*p[j]; // force evaluation
d=FABS(d);
if(d<delta)
d=delta;
tmp=p[j];
p[j]-=d;
(*func)(p, hxm, m, n, adata);
p[j]=tmp+d;
(*func)(p, hxp, m, n, adata);
p[j]=tmp; /* restore */
d=CNST(0.5)/d; /* invert so that divisions can be carried out faster as multiplications */
for(i=0; i<n; ++i){
jac[i*m+j]=(hxp[i]-hxm[i])*d;
}
}
}
static double evaluateCAT_FLEX (int *cptr, int *wptr,
double *x1, double *x2, double *tipVector,
unsigned char *tipX1, int n, double *diagptable_start, const int states)
{
double
sum = 0.0,
term,
*diagptable,
*left,
*right;
int
i,
l;
/* chosing between tip vectors and non tip vectors is identical in all flavors of this function ,regardless
of whether we are using CAT, GAMMA, DNA or protein data etc */
if(tipX1)
{
for (i = 0; i < n; i++)
{
/* same as in the GAMMA implementation */
left = &(tipVector[states * tipX1[i]]);
right = &(x2[states * i]);
/* important difference here, we do not have, as for GAMMA
4 P matrices assigned to each site, but just one. However those
P-Matrices can be different for the sites.
Hence we index into the precalculated P-matrices for individual sites
via the category pointer cptr[i]
*/
diagptable = &diagptable_start[states * cptr[i]];
/* similar to gamma, with the only difference that we do not integrate (sum)
over the discrete gamma rates, but simply compute the likelihood of the
site and the given P-matrix */
for(l = 0, term = 0.0; l < states; l++)
term += left[l] * right[l] * diagptable[l];
/* take the log */
term = LOG(FABS(term));
/*
multiply the log with the pattern weight of this site.
The site pattern for which we just computed the likelihood may
represent several alignment columns sites that have been compressed
into one site pattern if they are exactly identical AND evolve under the same model,
i.e., form part of the same partition.
*/
sum += wptr[i] * term;
}
}
else
{
for (i = 0; i < n; i++)
{
/* as before we now access the likelihood arrayes of two inner nodes */
left = &x1[states * i];
right = &x2[states * i];
diagptable = &diagptable_start[states * cptr[i]];
for(l = 0, term = 0.0; l < states; l++)
term += left[l] * right[l] * diagptable[l];
term = LOG(FABS(term));
sum += wptr[i] * term;
}
}
return sum;
}
/*
* This function returns the solution of Ax = b
*
* The function employs LU decomposition followed by forward/back substitution (see
* also the LAPACK-based LU solver above)
*
* A is mxm, b is mx1
*
* The function returns 0 in case of error, 1 if successful
*
* This function is often called repetitively to solve problems of identical
* dimensions. To avoid repetitive malloc's and free's, allocated memory is
* retained between calls and free'd-malloc'ed when not of the appropriate size.
* A call with NULL as the first argument forces this memory to be released.
*/
int AX_EQ_B_LU(LM_REAL *A, LM_REAL *B, LM_REAL *x, int m)
{
__STATIC__ void *buf=NULL;
__STATIC__ int buf_sz=0;
register int i, j, k;
int *idx, maxi=-1, idx_sz, a_sz, work_sz, tot_sz;
LM_REAL *a, *work, max, sum, tmp;
if(!A)
#ifdef LINSOLVERS_RETAIN_MEMORY
{
if(buf) free(buf);
buf=NULL;
buf_sz=0;
return 1;
}
#else
return 1; /* NOP */
#endif /* LINSOLVERS_RETAIN_MEMORY */
/* calculate required memory size */
idx_sz=m;
a_sz=m*m;
work_sz=m;
tot_sz=(a_sz+work_sz)*sizeof(LM_REAL) + idx_sz*sizeof(int); /* should be arranged in that order for proper doubles alignment */
// Check inputs for validity
for(i=0; i<a_sz; i++)
if (!LM_FINITE(A[i]))
return 0;
#ifdef LINSOLVERS_RETAIN_MEMORY
if(tot_sz>buf_sz){ /* insufficient memory, allocate a "big" memory chunk at once */
if(buf) free(buf); /* free previously allocated memory */
buf_sz=tot_sz;
buf=(void *)malloc(tot_sz);
if(!buf){
fprintf(stderr, RCAT("memory allocation in ", AX_EQ_B_LU) "() failed!\n");
exit(1);
}
}
#else
buf_sz=tot_sz;
buf=(void *)malloc(tot_sz);
if(!buf){
fprintf(stderr, RCAT("memory allocation in ", AX_EQ_B_LU) "() failed!\n");
exit(1);
}
#endif /* LINSOLVERS_RETAIN_MEMORY */
a=(LM_REAL*)buf;
work=a+a_sz;
idx=(int *)(work+work_sz);
/* avoid destroying A, B by copying them to a, x resp. */
memcpy(a, A, a_sz*sizeof(LM_REAL));
memcpy(x, B, m*sizeof(LM_REAL));
/* compute the LU decomposition of a row permutation of matrix a; the permutation itself is saved in idx[] */
for(i=0; i<m; ++i){
max=0.0;
for(j=0; j<m; ++j)
if((tmp=FABS(a[i*m+j]))>max)
max=tmp;
if(max==0.0){
fprintf(stderr, RCAT("Singular matrix A in ", AX_EQ_B_LU) "()!\n");
#ifndef LINSOLVERS_RETAIN_MEMORY
free(buf);
#endif
return 0;
}
work[i]=LM_CNST(1.0)/max;
}
for(j=0; j<m; ++j){
for(i=0; i<j; ++i){
sum=a[i*m+j];
for(k=0; k<i; ++k)
sum-=a[i*m+k]*a[k*m+j];
a[i*m+j]=sum;
//.........这里部分代码省略.........
static void
LNSRCH(int m, LM_REAL *x, LM_REAL f, LM_REAL *g, LM_REAL *p, LM_REAL alpha, LM_REAL *xpls,
LM_REAL *ffpls, void (*func)(LM_REAL *p, LM_REAL *hx, int m, int n, void *adata), struct FUNC_STATE state,
int *mxtake, int *iretcd, LM_REAL stepmx, LM_REAL steptl, LM_REAL *sx)
{
/* Find a next newton iterate by backtracking line search.
* Specifically, finds a \lambda such that for a fixed alpha<0.5 (usually 1e-4),
* f(x + \lambda*p) <= f(x) + alpha * \lambda * g^T*p
*
* Translated (with minor changes) from Schnabel, Koontz & Weiss uncmin.f, v1.3
* PARAMETERS :
* m --> dimension of problem (i.e. number of variables)
* x(m) --> old iterate: x[k-1]
* f --> function value at old iterate, f(x)
* g(m) --> gradient at old iterate, g(x), or approximate
* p(m) --> non-zero newton step
* alpha --> fixed constant < 0.5 for line search (see above)
* xpls(m) <-- new iterate x[k]
* ffpls <-- function value at new iterate, f(xpls)
* func --> name of subroutine to evaluate function
* state <--> information other than x and m that func requires.
* state is not modified in xlnsrch (but can be modified by func).
* iretcd <-- return code
* mxtake <-- boolean flag indicating step of maximum length used
* stepmx --> maximum allowable step size
* steptl --> relative step size at which successive iterates
* considered close enough to terminate algorithm
* sx(m) --> diagonal scaling matrix for x, can be NULL
* internal variables
* sln newton length
* rln relative length of newton step
*/
register int i;
int firstback = 1;
LM_REAL disc;
LM_REAL a3, b;
LM_REAL t1, t2, t3, lambda, tlmbda, rmnlmb;
LM_REAL scl, rln, sln, slp;
LM_REAL tmp1, tmp2;
LM_REAL fpls, pfpls = 0., plmbda = 0.; /* -Wall */
f*=CNST(0.5);
*mxtake = 0;
*iretcd = 2;
tmp1 = 0.;
if(!sx) /* no scaling */
for (i = 0; i < m; ++i)
tmp1 += p[i] * p[i];
else
for (i = 0; i < m; ++i)
tmp1 += sx[i] * sx[i] * p[i] * p[i];
sln = (LM_REAL)sqrt(tmp1);
if (sln > stepmx) {
/* newton step longer than maximum allowed */
scl = stepmx / sln;
for(i=0; i<m; ++i) /* p * scl */
p[i]*=scl;
sln = stepmx;
}
for(i=0, slp=0.; i<m; ++i) /* g^T * p */
slp+=g[i]*p[i];
rln = 0.;
if(!sx) /* no scaling */
for (i = 0; i < m; ++i) {
tmp1 = (FABS(x[i])>=CNST(1.))? FABS(x[i]) : CNST(1.);
tmp2 = FABS(p[i])/tmp1;
if(rln < tmp2) rln = tmp2;
}
else
for (i = 0; i < m; ++i) {
tmp1 = (FABS(x[i])>=CNST(1.)/sx[i])? FABS(x[i]) : CNST(1.)/sx[i];
tmp2 = FABS(p[i])/tmp1;
if(rln < tmp2) rln = tmp2;
}
rmnlmb = steptl / rln;
lambda = CNST(1.0);
/* check if new iterate satisfactory. generate new lambda if necessary. */
while(*iretcd > 1) {
for (i = 0; i < m; ++i)
xpls[i] = x[i] + lambda * p[i];
/* evaluate function at new point */
(*func)(xpls, state.hx, m, state.n, state.adata);
for(i=0, tmp1=0.0; i<state.n; ++i){
state.hx[i]=tmp2=state.x[i]-state.hx[i];
tmp1+=tmp2*tmp2;
}
fpls=CNST(0.5)*tmp1; *ffpls=tmp1; ++(*(state.nfev));
if (fpls <= f + slp * alpha * lambda) { /* solution found */
*iretcd = 0;
if (lambda == CNST(1.) && sln > stepmx * CNST(.99)) *mxtake = 1;
return;
//.........这里部分代码省略.........
客服电话
APP下载
官方微信
请发表评论