[xiph-commits] r17898 - trunk/ghost/monty/lift
xiphmont at svn.xiph.org
xiphmont at svn.xiph.org
Fri Mar 18 00:09:21 PDT 2011
Author: xiphmont
Date: 2011-03-18 00:09:21 -0700 (Fri, 18 Mar 2011)
New Revision: 17898
Added:
trunk/ghost/monty/lift/anneal.c
Modified:
trunk/ghost/monty/lift/granule_simplex.c
trunk/ghost/monty/lift/liftsearch.c
Log:
Code cleanup / code catchup
Added: trunk/ghost/monty/lift/anneal.c
===================================================================
--- trunk/ghost/monty/lift/anneal.c (rev 0)
+++ trunk/ghost/monty/lift/anneal.c 2011-03-18 07:09:21 UTC (rev 17898)
@@ -0,0 +1,159 @@
+/********************************************************************
+ * *
+ * THIS FILE IS PART OF THE OggGhost SOFTWARE CODEC SOURCE CODE. *
+ * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS *
+ * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
+ * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. *
+ * *
+ * THE OggGhost SOURCE CODE IS (C) COPYRIGHT 1994-2007 *
+ * the Xiph.Org FOundation http://www.xiph.org/ *
+ * *
+ ********************************************************************
+
+ function: implement adaptive simulated annealing
+ last mod: $Id$
+
+ ********************************************************************/
+
+/* see: _Digital IIR Filter Design Using Adaptive Simpulated
+ Annealing_ by S. Chen, R. Istepanian, and B. L. Luk */
+
+#include <stdlib.h>
+#include <unistd.h>
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+
+static double g(double Ti){
+ double ui = drand48();
+ if(ui<.5)
+ return -Ti * (pow(1.+1./Ti, fabs(2*ui-1))-1);
+ else
+ return Ti * (pow(1.+1./Ti, fabs(2*ui-1))-1);
+}
+
+static inline double todB(double x){
+ return ((x)==0?-400:log((x)*(x))*4.34294480f);
+}
+
+void anneal(double *w, int m,
+ double c, double U, double V, double epsilon,
+ int N_accept, int N_generate, int N_break,
+ double (*cost_func)(double *)){
+
+ // (i)
+
+ double cost = cost_func(w);
+ double cost_best = cost;
+ double w_best[m];
+ double Tc0 = cost;
+ double Tc = Tc0;
+ double kc = 0.;
+ double T[m];
+ double k[m];
+ int accepted = 0;
+ int generated = 0;
+ int bests = 0;
+ int last_bests = 0;
+ int break_counter = 0;
+ int i;
+
+ memcpy(w_best,w,sizeof(w_best));
+ for(i=0;i<m;i++){
+ T[i] = 1.;
+ k[i] = 0.;
+ }
+
+ fprintf(stderr,"\nannealing... %d/%d break:%d cost:%g ",
+ accepted,generated,N_break-break_counter,todB(cost_best)/2);
+
+ while(break_counter<N_break){
+
+ // (ii) generate a new point
+
+ double w_new[m];
+ for(i=0; i<m; i++){
+ while(1){
+ w_new[i] = w[i] + g(T[i])*(V-U);
+ if(w_new[i]<U)continue;
+ if(w_new[i]>V)continue;
+ break;
+ }
+ }
+
+ double cost_new = cost_func(w_new);
+ double P_accept = 1. / (1. + exp( (cost_new - cost)/Tc));
+ double P_unif = drand48();
+
+ if(P_unif <= P_accept){
+ cost = cost_new;
+ memcpy(w,w_new,sizeof(w_new));
+ accepted++;
+
+ if(cost < cost_best){
+ cost_best = cost;
+ memcpy(w_best,w,sizeof(w_best));
+ bests++;
+
+ fprintf(stderr,"\rannealing... %d/%d break:%d cost:%g ",
+ accepted,generated,N_break-break_counter,todB(cost_best)/2);
+ }
+ }
+ generated++;
+
+ // (iii) reannealing
+ if(accepted >= N_accept){
+ double s[m];
+ for(i=0; i<m; i++){
+ memcpy(w_new,w_best,sizeof(w_best));
+ w_new[i]+=epsilon;
+
+ s[i] = fabs((cost_func(w_new) - cost_best) / epsilon);
+ }
+ double s_max = s[0];
+ for(i=1; i<m; i++)
+ if(s_max < s[i]) s_max = s[i];
+
+ for(i=0;i<m;i++){
+ T[i] *= s_max/s[i];
+ k[i] = pow(-log(T[i])/c, m);
+ }
+
+ Tc0 = cost;
+ Tc = cost_best;
+ kc = pow(-log(Tc/Tc0)/c, m);
+
+ if(last_bests == bests)
+ break_counter++;
+ else
+ break_counter=0;
+
+ last_bests = bests;
+
+ fprintf(stderr,"\rannealing... %d/%d break:%d cost:%g ",
+ accepted,generated,N_break-break_counter,todB(cost_best)/2);
+ accepted = 0;
+ }
+
+ // (iv) temperature annealing
+ if(generated >= N_generate){
+
+ for(i=0;i<m;i++){
+ k[i] += 1.;
+ T[i] = exp(-c * pow(k[i], 1./m));
+ }
+
+ kc += 1.;
+ Tc = Tc0 * exp(-c * pow(kc, 1./m));
+
+ fprintf(stderr,"\rannealing... %d/%d break:%d cost:%g ",
+ accepted,generated,N_break-break_counter,todB(cost_best)/2);
+ generated=0;
+ }
+ }
+
+ fprintf(stderr,"\rannealed. %d/%d, cost:%g ",
+ accepted,generated,todB(cost_best)/2);
+
+ memcpy(w,w_best,sizeof(w_best));
+}
Modified: trunk/ghost/monty/lift/granule_simplex.c
===================================================================
--- trunk/ghost/monty/lift/granule_simplex.c 2011-03-18 06:42:25 UTC (rev 17897)
+++ trunk/ghost/monty/lift/granule_simplex.c 2011-03-18 07:09:21 UTC (rev 17898)
@@ -99,19 +99,148 @@
return cost;
}
-void grow_vertex(int *c, int m,
- int *allowed_movements,
+int grow_vertex(int *c, int m,
+ int *movement,
int min_bound, int max_bound,
double (*cost_func)(int *)){
+
+ int changeflag = 1;
+ int i,count=0;
+ double cost = cost_func(c);
+ while(changeflag){
+ changeflag=0;
+ for(i=0;i<m;i++){
+ if(movement[i]){
+ c[i] += movement[i];
+ double test_cost = lift_cost(c);
+ if(test_cost>cost){
+ changeflag=1;
+ cost=test_cost;
+ count++;
+ }else{
+ c[i] -= movement[i];
+ }
+ }
+ }
+ }
+ return count;
+}
+void center_of_tope(int *center, int **verticies, int m){
+ int acc[m];
+ int i, j;
+ memset(acc,0,sizeof(acc));
+
+ for(i=0;i<m-1;i++)
+ for(j=0;j<m;j++)
+ acc[j]+=verticies[i][j];
+
+ for(j=0;j<m;j++)
+ center[j] = acc[j] / (m-1);
}
+void grow_secondary_granule(int *minimumc,
+ int **c, int m,
+ int *movemask,
+ int min_bound, int max_bound,
+ double (*cost_func)(int *),
+ int *simplexen,
+ int *granules){
+ int i;
+
+ // grow the single free vertex
+ int movement = grow_vertex(c[0], m, movemask, min_bound, max_bound, cost_func);
+
+ if(movement<=1){ // 1 is roundoff error slack
+ // we cannot grow a new simplex; the simplex facet is pressed flat
+ // against a convergence ridge.
+
+ // this is the correct point to crest the ridge and walk the granule
+ // down the other side looking for a new minimum. If it is indeed a new
+ // minimum, grow a new granule.
+
+ for(i=0;i<m;i++)
+ c[0][i] += movemask[i];
+
+ double cost = gradient_walk(c[0], m, min_bound, max_bound, cost_func);
+ int min_num = lookup_minimum(c[0]);
+ if(min_num!=-1){
+ log_minimum(c[0],cost);
+ grow_new_granule(c[0]); // enter recursion
+ }
+
+ return;
+ }
+
+ // continue recursively spawning secondary simplexes with one
+ // unconstrained vertex from the center of each face of this
+ // simplex. The new vertex must grow outward (reflect the simplex
+ // through each face).
+ int **secondary_c = malloc((m+1) * sizeof(*secondary_c));
+ secondary_c[0] = malloc(m * sizeof(**new_c));
+
+ for(i=0;i<(m+1);i++){
+ // build a new d-1-tope 'plane'; all the points in the current
+ // simplex minus the i'th one.
+
+ // fill in the fixed verticies [1 through m]
+ for(j=1;j<m;j++)
+ if(j<=i)
+ secondary_c[j] = new_c[j-1];
+ else
+ secondary_c[j] = new_c[j];
+
+ // first vertex is the unconstrained one
+ center_of_tope(secondary_c[0], secondary_c+1, m);
+
+ // set movemask [reverse of the i'th vertex's movemask]
+ memset(movemask,0,sizeof(*movemask)*m);
+ if(i>0)movemask[i-1]=1;
+ if(i<m)movemask[i]=-1;
+
+ grow_secondary_granule(secondary_c, m, movemask, min_bound, max_bound, cost_func, simplexen);
+ }
+
+ // Each vertex is sitting on a ridge. Pop over the ridge and
+ // walk gradient down to a minimum. If it's a new minimum, grow a
+ // new granule.
+
+ // this is here and not above to allow each granule to fill out
+ // completely before beginning work on another.
+
+
+ for(i=0;i<(m+1);i++){
+ // re-set up orthogonal move mask for this vertex
+ memset(movemask,0,sizeof(*movemask)*m);
+ if(i>0)movemask[i-1]=-1;
+ if(i<m)movemask[i]=1;
+ //crest_walk_and_recurse(new_c[i], m, movemask, min_bound, max_bound, cost_func, granules);
+ free(new_c[i]);
+ }
+
+ free(movemask);
+ free(new_c);
+
+
+
+
+
+
+
+
+
+}
+
+
+
void grow_new_granule(int *c, int m,
int min_bound, int max_bound,
double (*cost_func)(int *)){
int i;
+ int simplexen = 1;
+ int granules = 1;
// Start by growing a new simplex by allowing each vertex to move
// outward from the center (the local minimum). A vertex continues
// moving so long as it is moving uphill. For this new simplex, the
@@ -155,12 +284,11 @@
// first vertex is the unconstrained one
center_of_tope(secondary_c[0], secondary_c+1, m);
- // set movemask [reverse of the i'th vertex's movemask]
- memset(movemask,0,sizeof(*movemask)*m);
- if(i>0)movemask[i-1]=1;
- if(i<m)movemask[i]=-1;
+ // set movemask; must be strictly away from local minimum
+ movemask_away(c, movemask, m);
- grow_secondary_granule(secondary_c, m, movemask, min_bound, max_bound, cost_func);
+ // grow granule
+ grow_secondary_granule(secondary_c, m, movemask, min_bound, max_bound, cost_func, simplexen);
}
// Each vertex is sitting on a ridge. Pop over the ridge and
@@ -169,12 +297,14 @@
// this is here and not above to allow each granule to fill out
// completely before beginning work on another.
+
+
for(i=0;i<(m+1);i++){
// re-set up orthogonal move mask for this vertex
memset(movemask,0,sizeof(*movemask)*m);
if(i>0)movemask[i-1]=-1;
if(i<m)movemask[i]=1;
- crest_walk_and_recurse(new_c[i], m, movemask, min_bound, max_bound, cost_func);
+ //crest_walk_and_recurse(new_c[i], m, movemask, min_bound, max_bound, cost_func, granules);
free(new_c[i]);
}
Modified: trunk/ghost/monty/lift/liftsearch.c
===================================================================
--- trunk/ghost/monty/lift/liftsearch.c 2011-03-18 06:42:25 UTC (rev 17897)
+++ trunk/ghost/monty/lift/liftsearch.c 2011-03-18 07:09:21 UTC (rev 17898)
@@ -441,7 +441,7 @@
void set_maxflat(double *fc){
int i, j;
for(i=0,j=0;i<order;i++,j++)
- fc[j] = def_a(order, i, order-1)+(drand48()-.5)*.01;
+ fc[j] = def_a(order, i, order-1);//+(drand48()-.5)*.01;
//for(i=0;i<order;i++,j++)
// fc[j] = def_a(order, order-i-1, order-1)+(drand48()-.5)*.01;
@@ -453,6 +453,7 @@
//fc[j++]=1.;
}
+
int main(int argc, char *argv[]){
drft_init(&fft,FS);
@@ -466,10 +467,12 @@
s_delay = 0;
//double *fc = calloc(num_coefficients, sizeof(*fc));
+ //set_maxflat(fc);
- //set_maxflat(fc);
+<<<<<<< .mine
+ //lift_dump(fc);
+ //exit(0);
//num_coefficients = (order*2+1)*2;
-
//walk_to_minimum_A(fc,1./(1<<10));
//walk_to_minimum_A(fc,1./(1<<12));
//walk_to_minimum_A(fc,1./(1<<14));
@@ -478,15 +481,25 @@
//lift_dump(fc);
//fflush(stdout);
+ //double fc[16]={ 0.494814, -0.117399, 0.053387, -0.029044, 0.016904, -0.010041,
+ // 0.005928, -0.003415, 0.001890, -0.000988, 0.000478, -0.000208,
+ // 0.000078, -0.000022, 0.000003, 0.000001 };
- double fc[16]={ 0.494814, -0.117399, 0.053387, -0.029044, 0.016904, -0.010041,
- 0.005928, -0.003415, 0.001890, -0.000988, 0.000478, -0.000208,
- 0.000078, -0.000022, 0.000003, 0.000001 };
+ //double fc[16] = {0.496302, -0.119542, 0.055873, -0.031647, 0.019441, -0.012371, 0.007955, -0.005087,
+ // 0.003196, -0.001951, 0.001144, -0.000635, 0.000327, -0.000152, 0.000060, -0.000017};
+ double fc[16] = {0.496355, -0.1196186, 0.05596481, -0.03174554, 0.01954113, -0.01246747, 0.008044547, -0.005165739,
+ 0.003263023, -0.002006148, 0.001186254, -0.0006671808, 0.0003496294, -0.0001656954, 0.00006832193, -0.00002147376};
+
+ //double fc[16]= {0.4970733, -0.1206598, 0.05718326, -0.03303792,
+ // 0.02083225, -0.0136748, 0.009122373, -0.006084431,
+ // 0.004005327, -0.002563749, 0.001590537, -0.0009246328,
+ // 0.0005063734, -0.000228745, 8.170831e-05, -2.94743e-05};
+
//set_maxflat(fc);
//walk_to_minimum_CSD(fc,1./(1<<16));
- //lift_dump(fc);
- //fflush(stdout);
+ lift_dump(fc);
+ fflush(stdout);
//walk_to_minimum_CSD(fc,1./(1<<17));
//lift_dump(fc);
//fflush(stdout);
@@ -513,7 +526,6 @@
//lift_dump(fc);
//fflush(stdout);
-
lift_dump(fc);
return 0;
More information about the commits
mailing list