http://www.ingber.com/asa_contrib.txt ftp://ftp.ingber.com/asa_contrib.txt This file contains code contributed or developed by other authors that is directly useful for ASA. Some "toy" problems optimized using ASA, which can provide immediate examples on how you can optimize your own problem, are given in http://www.ingber.com/asa_examples.txt ftp://ftp.ingber.com/asa_examples.txt For papers dealing with difficult systems across a wide range of disciplines, see http://www.ingber.com/asa_papers.html ftp://ftp.ingber.com/asa_papers.txt CONTENTS @@CHANGES @@ASA DLL @@ASA_PARALLEL @@MATLAB Interface with ASA @@RLAB Interface with ASA @@AMPL Interface to ASA @@Haskell Interface with ASA @@FUZZY_ASA @@CHIPSET for ASA_QUEUE @@Hashed Lists for ASA_QUEUE @@FDLIBM Code for FDLIBM_EXP, FDLIBM_LOG and FDLIBM_POW @@RATFOR vfsr.r Code (1987) unsupported ======================================================================== ======================================================================== @@CHANGES 22 Feb 03 Renamed files to faciliate their use in larger projects: asa_user.h -> asa_usr_asa.h user.c -> asa_usr.c user.h -> asa_usr.h user_cst.c -> asa_usr_cst.c readme.ms -> ASA-README.ms test_asa -> asa_test_asa test_usr -> asa_test_usr ======================================================================== @@ASA DLL 10 Apr 03 Creating an ASA DLL is quite easy under Cygwin (cygwin.com) on Windows. Copy the ASA distributiuon into a new directory, e.g., ASA_DLL. (1) Modify ASA with this simple patch (as of Version 25.1) ------------8<------------ top cut -> bottom ------------->8------------ diff -rc ASA/asa_usr.c /usr/local/doc/DLL/ASA_DLL/asa_usr.c *** ASA/asa_usr.c Sat Feb 22 07:34:48 2003 --- /usr/local/doc/DLL/ASA_DLL/asa_usr.c Wed Apr 9 10:54:57 2003 *************** *** 4628,4633 **** --- 4628,4634 ---- } #endif /* ASA_TEMPLATE_SAMPLE */ + #if 0 #if ASA_TEMPLATE_LIB int main () *************** *** 4675,4680 **** --- 4676,4682 ---- /* NOTREACHED */ } #endif /* ASA_TEMPLATE_LIB */ + #endif /* if 0 */ void Exit_USER (char *statement) Only in /usr/local/doc/DLL/ASA_DLL/: asa_usr.o diff -rc ASA/asa_usr_asa.h /usr/local/doc/DLL/ASA_DLL/asa_usr_asa.h *** ASA/asa_usr_asa.h Sat Feb 22 07:34:55 2003 --- /usr/local/doc/DLL/ASA_DLL/asa_usr_asa.h Wed Apr 9 11:01:08 2003 *************** *** 41,46 **** --- 41,47 ---- #endif #if MY_TEMPLATE /* MY_TEMPLATE_asa_user */ /* you can add your own set of #define here */ + #define ASA_TEMPLATE_LIB TRUE #endif /* MY_TEMPLATE */ #ifndef ASA_TEMPLATE_LIB ------------8<------------ bottom cut <- top ------------->8------------ (2) Create a main.c in your new directory, essentially the main() in asa_usr.c used for ASA_TEMPLATE_LIB: ------------8<------------ top cut -> bottom ------------->8------------ #include "asa_usr.h" #if ASA_TEMPLATE_LIB int main () { double main_cost_value; double *main_cost_parameters; int main_exit_code; LONG_INT number_params; ALLOC_INT n_param; FILE *ptr_main; char user_exit_msg[160]; ptr_main = stdout; number_params = 4; if ((main_cost_parameters = (double *) calloc (number_params, sizeof (double))) == NULL) { strcpy (user_exit_msg, "ASA_TEMPLATE_LIB main(): main_cost_parameters"); Exit_USER (user_exit_msg); return (-2); } asa_seed (696969); /* This is the default random seed. */ asa_main (&main_cost_value, main_cost_parameters, &main_exit_code); fprintf (ptr_main, "main_exit_code = %d\n", main_exit_code); fprintf (ptr_main, "main_cost_value = %12.7g\n", main_cost_value); fprintf (ptr_main, "parameter\tvalue\n"); for (n_param = 0; n_param < number_params; ++n_param) { fprintf (ptr_main, "%d\t\t%12.7g\n", n_param, main_cost_parameters[n_param]); } free (main_cost_parameters); return (0); } #endif /* ASA_TEMPLATE_LIB */ ------------8<------------ bottom cut <- top ------------->8------------ (3) Run this script: ------------8<------------ top cut -> bottom ------------->8------------ #!/bin/tcsh -fxv set Cfiles = ( asa asa_usr asa_usr_cst ) foreach i ( $Cfiles ) gcc -mrtd -mno-cygwin -o $i.o -c $i.c end set module = asadll set old_lib = " asa.o asa_usr.o asa_usr_cst.o " set dependency_libs = '' gcc -mrtd -shared -mno-cygwin -o ${module}.dll \ -Wl,--out-implib=lib${module}.a \ -Wl,--export-all-symbols \ -Wl,--export-dynamic \ -Wl,--add-stdcall-alias \ -Wl,--enable-auto-import \ -Wl,--output-def=${module}.def \ -Wl,--whole-archive ${old_lib} \ -Wl,--no-whole-archive ${dependency_libs} gcc -mrtd -mno-cygwin -O -o run_main main.c -L./ -l${module} ${dependency_libs} ./run_main ------------8<------------ bottom cut <- top ------------->8------------ Lester ======================================================================== @@ASA_PARALLEL 24 Jan 01 19 Sep 06 In the early 1990's I used the the ASA_PARALLEL hooks with success on a Connection machine and in 1994-1995 on Crays for an NSF project, Parallelizing ASA and PATHINT Project (PAPP). In 2001 Thomas Binder developed auxiliary functions to support ASA_PARALLEL for ASA version 1.15. ======================================================================== @@MATLAB Interface with ASA 6 Jan 03 Note the section Generic ASA Interfaces in ASA-README, which outlines how to generically interface ASA to other languages. I have used this as well for a MATLAB interface. The contribution below will be useful for a standalone simpler approach to using ASA with MATLAB, albeit with some limitations on the use of many ASA OPTIONS. 8 Dec 99 [address updated 12 Aug 03] To: ingber@ingber.com Subject: MATLAB gateway routine to ASA From: Shinichi Sakata Dear Lester, Asamin provides an interface (mex file) for using ASA within MATLAB. With asamin, you call ASA from MATLAB to search the minimum of an objective function coded in the MATLAB language. Asamin also offers a way to change many of ASA's run-time options. You would not need much knowledge about the C language to use asamin, unless you want to use an array of compilation-time options available for ASA. Asamin is written for and tested with MATLAB 5.3 on Solaris 2.5.1. In order to run it with other versions of MATLAB or on other platforms, you may need to make some changes to its source code. Asamin is available from 1 Jun 04 [updated] http://www.econ.ubc.ca/ssakata/public_html/software/ Have a happy holiday, Shinichi 22 Nov 07 You should be able to control most ASA post-compile and pre-compile options by setting them in asa_usr_asa.h. If some of those you require are set in asa_usr.c, and not in asamin.c (which replaces asa_usr.c and asa_usr.h), you would have to add some variables to asamin.c and asamin.h (or see alternative below). You do not have to be concerned with platforms that must read in a parameter file like asa_opt. In other words, some options that are in the standard ASA-Makefile can instead be set in asa_usr_asa.h (see http://www.ingber.com/asa_examples.txt), and if OPTIONS_FILE is set to FALSE there as well, then all the options in asa_opt (and more) can be set directly in asamin.c, as they might appear in asa_usr.c (where there are ASA_TEMPLATEs provided to illustrate how to set them). Alternatively, include asa_usr.c and asa_usr.h, and set up a small mex file similar to asamin.c, but which calls asa_main() in asa_usr.c, and pass only those variables between Matlab <-> asamin() <-> asa_main() as required. You would have to re-compile these C codes if you change any of the "hard-coded" OPTIONS in the ASA files. So, if you are using the Matlab asamin interface, AND you really truly honestly believe that your project really truly honestly requires some options not "obviously" already available in asamin, you can try the above advanced procedures. Also note the section "Use of Documentation for Tuning" in ASA-README-.[] which may be vital for use of ASA (or in fact for just about any sampling algorithm powerful enough to apply to many classes of hard problems -- since such problems are typically non-typical and some tuning hard work often is required). Lester ======================================================================== @@RLAB Interface with ASA 18 Feb 08 Note that RLaB's asamin() is not the same function as the Matlab asamin() discussed above. They both have developed a user module using the ASA asa module. Lester Date: Mon, 18 Feb 2008 14:43:44 -0500 From: mkostrun@gmail.com To: ingber@ingber.com Subject: asa code Dear Dr. Ingber, this is to inform you that I have incorporated your ASA code into the project rlabplus (matlab-like environment), http://rlabplus.sourceforge.net, where it is now accessible as a built-in function asamin(..). The source code distribution of the rlabplus contains the following files asa.c, asa.h, usr_asa_usr.h and asa-license.txt, neither of which was modified. The file Makefile.in, which is a build file for rlabplus, contains a list of define's particular to my implementation. I am grateful that you have made your excellent solver available to use in public domain. Regards, Marijan Kostrun ============================================================================== @@AMPL Interface to ASA 22 Apr 07 Ismael Vaz has written an AMPL interface to ASA. This allows the ASA solver to be available through the NEOS server at http://neos.mcs.anl.gov/neos/solvers/go:ASA/AMPL.html. ============================================================================== @@Haskell Interface with ASA 25 Jan 06 John Meacham has written a Haskell [http://haskell.org/] interface with ASA, HsASA. His Docs and source code can be downloaded from http://repetae.net/john/recent/out/HsASA.html . ============================================================================== @@FUZZY_ASA 28 Jun 04 The link the the paper below is now http://www.optimization-online.org/DB_FILE/2003/11/779.pdf 6 Jan 03 Note the section Generic ASA Interfaces in ASA-README, which outlines how to generically interface ASA to other languages. I have used this as well for a interfaces to complex cost functions. The contribution below will be useful for a standalone simple approach to using ASA with a specific algorithm to adaptively change an ASA OPTIONS, albeit with some limitations on the use of many ASA OPTIONS. 29 Oct 02 Here is an interesting application of a fuzzy algorithm to adaptively tune a given parameter set of ASA OPTIONS. This approach has been developed by %A H.A.O. Junior %T Fuzzy control of stochastic global optimization algorithms and very fast simulated reannealing %I hime@engineer.com %C Rio de Janeiro, Brazil %D 2001 Note that MULTI_MIN also might be used to generate a set of points to which this kind of fuzzy algorithm can be applied. The following three files, ExampleCost.c, FuzzyASAExample.c, and fuzzyasaexample.bat were contributed by Hime. ------------8<------------ top cut -> bottom ------------->8------------ /* :::::::::::::: ExampleCost.c :::::::::::::: */ #include ///////// This example function was extracted from the book ///////// Global optimization in Action by J nos D. Pint‚r ///////// Kluwer Academic Publishers - 1996 double ExampleCost(double *x,int nodedimensoes) { int i; double f1,f2,f3,aux; f1=0; for (i=0;iLimit_Acceptances = 1000000; USER_OPTIONS->Limit_Generated = 2000000; USER_OPTIONS->Limit_Invalid_Generated_States = 100000000; USER_OPTIONS->Accepted_To_Generated_Ratio = 1.0E-4; USER_OPTIONS->Cost_Precision = 1.0E-18; USER_OPTIONS->Maximum_Cost_Repeat = 350; // ERA 35 USER_OPTIONS->Number_Cost_Samples = 25; USER_OPTIONS->Temperature_Ratio_Scale = .9 ; // Very important parameter ! USER_OPTIONS->Cost_Parameter_Scale_Ratio = 1.0; USER_OPTIONS->Temperature_Anneal_Scale = 100.0; USER_OPTIONS->Include_Integer_Parameters = FALSE; USER_OPTIONS->User_Initial_Parameters = TRUE; USER_OPTIONS->Sequential_Parameters = -1; USER_OPTIONS->Initial_Parameter_Temperature = 1000.0; USER_OPTIONS->Acceptance_Frequency_Modulus = 100; USER_OPTIONS->Generated_Frequency_Modulus = 10000; USER_OPTIONS->Reanneal_Cost = 1; USER_OPTIONS->Reanneal_Parameters = TRUE; USER_OPTIONS->Delta_X = 0.001; USER_OPTIONS->User_Tangents = FALSE; USER_OPTIONS->Curvature_0 = TRUE; USER_OPTIONS->User_Quench_Param_Scale = CoolingScales ; USER_OPTIONS->User_Quench_Cost_Scale = &ScaleOfCost ; for (index = 0 ; index < NOOFSAMPLES ; ++index) Parameters[index] = (double *) calloc(*parameter_dimension,sizeof(double)); cost_value = asa (MyCost, randflt, rand_seed, cost_parameters, parameter_lower_bound, parameter_upper_bound, cost_tangents, cost_curvature, parameter_dimension, parameter_int_real, cost_flag, exit_code, USER_OPTIONS); printf("\n Final value = %f \n",cost_value); for (index=0;index<*parameter_dimension;index++) printf("\n Parameter %d = %f",index+1,cost_parameters[index]); printf("\nTotal number of evaluations = %ld",Evaluations); } double MyCost (double *x, double *parameter_lower_bound, double *parameter_upper_bound, double *cost_tangents, double *cost_curvature, ALLOC_INT * parameter_dimension, int *parameter_int_real, int *cost_flag, int *exit_code, USER_DEFINES * USER_OPTIONS) { void Controla(double * , double , USER_DEFINES * , int ); double ExampleCost(double *x,int nodedimensoes); double valoratu=0; int i; Evaluations++; valoratu = ExampleCost(x,*parameter_dimension); *cost_flag = TRUE; //////////////////////////////////////////////////// switch( USER_OPTIONS->Locate_Cost ) { case -1: case 2: case 5: case 12: Controla(x,valoratu,USER_OPTIONS,*parameter_dimension); if (! (Evaluations % 100 ) ) fprintf(fp,"%ld %e\n",Evaluations,*(USER_OPTIONS->Best_Cost)); break; } //////////////////////////////////////////////////// return valoratu; } void MeanAndDeviation(double *Vetor , int NoOfElements , double *Media ,double *deviation ) { int i; double aux; *Media = *deviation = 0; for (i=0 ; i < NoOfElements ; i++) { aux = Vetor[i]; *Media += aux; *deviation += aux*aux; } *Media /= NoOfElements ; *deviation /= NoOfElements ; *deviation = *deviation - (*Media)*(*Media) ; *deviation = sqrt(fabs(*deviation)); } ///////// void Controla (double *x,double valorf,USER_DEFINES *USER_OPTIONS , int NoParam) { double SubEnergia(double Valor,double Minimo); void AlteraQuench(USER_DEFINES *opcoes,int NoParam,double Media,double Desvio); static double Desempenho , Media , Desvio , aux[NOOFSAMPLES] ; static int IndVal = 0; int i; IndVal++; if (IndVal % (NOOFSAMPLES+1)) { Values[IndVal] = valorf ; Minimums[IndVal] = * USER_OPTIONS->Best_Cost ; for (i=0;iBest_Parameters[i]; // Salva parametros return ; } printf("\nActivating controller"); IndVal=0; if (Minimums[1] !=0) Desempenho = ( Minimums[1] - Minimums[NOOFSAMPLES] )/fabs(Minimums[1]); else return; if ( Desempenho > META ) return ; for (i=0 ; i < NOOFSAMPLES ; i++) aux[i] = SubEnergia(Values[i+1],Minimums[i+1]); // Zero based MeanAndDeviation(aux,NOOFSAMPLES,&Media,&Desvio); AlteraQuench(USER_OPTIONS,NoParam,Media,Desvio); return; } double SubEnergia(double Valor,double Minimo) { #define a 2 double valuhat , argulog ; valuhat = Valor - Minimo ; argulog = 1/(1+exp( -(valuhat+a) ) ); if (argulog == 0) argulog = 1e-100; return log (argulog) ; } double FatorDelta(double MediaSub) { #define DELTATPOS .3 #define DELTATNULO 0 #define DELTATNEG -.3 #define VALMINLOC log(1/(1+exp(-a))) #define INTER1 .1 double PertiMediaZero =0 ; double PertiMediaMedia=0 ; double PertiMediaBaixa=0 ; double ResuCrisp; /////////////////// VERY SIMPLE FUZZY RULEBASE /////////////////// // Fuzzy rule #1 - IF MediaSub is ZERO THEN Decrease acceleration IS POSITIVE // ( We are in distant locations even from the present basic local minimum ) if (MediaSub > Threshold1) PertiMediaZero = (MediaSub - Threshold1)/(-Threshold1) ; // Fuzzy rule #2 - IF MediaSub IS VALMINLOC ( FUZZY NUMBER ) THEN // Decrease acceleration IS POSITIVE // ( We are in points near to present basic local minimum ) if ( MediaSub >= VALMINLOC && MediaSub <= VALMINLOC+INTER1 ) PertiMediaMedia = (VALMINLOC+INTER1 - MediaSub)/INTER1; ResuCrisp = PertiMediaZero*DELTATPOS + PertiMediaMedia*DELTATPOS ; return ResuCrisp; } void AlteraQuench(USER_DEFINES *opcoes , int NoParam , double Media,double deviation) { int i,j; double Delta,Mediaaux,deviationaux; double aux[NOOFSAMPLES]; Delta = FatorDelta(Media); if (opcoes->User_Quench_Cost_Scale[0] < 100) opcoes->User_Quench_Cost_Scale[0] *= (1 + .5*Delta) ; for (i=0;i < NoParam ;i++) if (opcoes->User_Quench_Param_Scale[i] < 100) opcoes->User_Quench_Param_Scale[i] *= (1 + .2*Delta) ; //0.05 if (deviation < THRESHOLDDEVIATION) { for (i=0;i < NoParam ;i++) if ( opcoes->User_Quench_Param_Scale[i] < 100 && opcoes->User_Quench_Param_Scale[i] > .1 ) { for (j=0;jUser_Quench_Param_Scale[i] /= (1+.15*Delta*exp(-deviationaux)); } } } //////////////////////////////////////////////////////////////////////////// /// The code below was extracted from USER.C (ASA distribution) //////////// //////////////////////////////////////////////////////////////////////////// /* Here is a good random number generator */ #define SHUFFLE 256 /* size of random array */ #define MULT ((LONG_INT) 25173) #define MOD ((LONG_INT) 65536) #define INCR ((LONG_INT) 13849) #define FMOD ((double) 65536.0) /*********************************************************************** * double myrand - returns random number between 0 and 1 * This routine returns the random number generator between 0 and 1 ***********************************************************************/ double myrand (LONG_INT * rand_seed) /* returns random number in {0,1} */ { #if TRUE /* (change to FALSE for alternative RNG) */ *rand_seed = (LONG_INT) ((MULT * (*rand_seed) + INCR) % MOD); return ((double) (*rand_seed) / FMOD); #else /* See "Random Number Generators: Good Ones Are Hard To Find," Park & Miller, CACM 31 (10) (October 1988) pp. 1192-1201. *********************************************************** THIS IMPLEMENTATION REQUIRES AT LEAST 32 BIT INTEGERS *********************************************************** */ #define _A_MULTIPLIER 16807L #define _M_MODULUS 2147483647L /* (2**31)-1 */ #define _Q_QUOTIENT 127773L /* 2147483647 / 16807 */ #define _R_REMAINDER 2836L /* 2147483647 % 16807 */ long lo; long hi; long test; hi = *rand_seed / _Q_QUOTIENT; lo = *rand_seed % _Q_QUOTIENT; test = _A_MULTIPLIER * lo - _R_REMAINDER * hi; if (test > 0) { *rand_seed = test; } else { *rand_seed = test + _M_MODULUS; } return ((double) *rand_seed / _M_MODULUS); #endif /* alternative RNG */ } /*********************************************************************** * double randflt ***********************************************************************/ double randflt (LONG_INT * rand_seed) /* shuffles random numbers in random_array[SHUFFLE] array */ { /* This RNG is a modified algorithm of that presented in * %A K. Binder * %A D. Stauffer * %T A simple introduction to Monte Carlo simulations and some * specialized topics * %B Applications of the Monte Carlo Method in statistical physics * %E K. Binder * %I Springer-Verlag * %C Berlin * %D 1985 * %P 1-36 * where it is stated that such algorithms have been found to be * quite satisfactory in many statistical physics applications. */ double rranf; unsigned kranf; int n; static int initial_flag = 0; static double random_array[SHUFFLE]; /* random variables */ if (initial_flag == 0 || *rand_seed < 0) { if (*rand_seed < 0) *rand_seed = -*rand_seed; for (n = 0; n < SHUFFLE; ++n) random_array[n] = myrand (rand_seed); initial_flag = 1; for (n = 0; n < 1000; ++n) /* warm up random generator */ rranf = randflt (rand_seed); return (rranf); } kranf = (unsigned) (myrand (rand_seed) * SHUFFLE) % SHUFFLE; rranf = *(random_array + kranf); *(random_array + kranf) = myrand (rand_seed); return (rranf); } /* :::::::::::::: fuzzyasaexample.bat :::::::::::::: */ rem Command line for compilation (Borland compiler). rem Adaptation to other compilers should be trivial. bcc32 -L\bcc55\lib -I\bcc55\include -w- -DASA_PRINT=FALSE -DNO_PARAM_TEMP_TEST=TRUE -DNO_COST_TEMP_TEST=TRUE -DQUENCH_COST=TRUE -DQUENCH_PARAMETERS=TRUE -DMAXIMUM_REANNEAL_INDEX=100000 -DCHECK_EXPONENT=TRUE -DMY_TEMPLATE=FALSE -DSMALL_FLOAT=1e-300 -DOPTIONS_FILE=FALSE -DOPTIONS_FILE_DATA=FALSE -DCOST_FILE=FALSE -DASA_LIB=TRUE -eFuzzyASAExample FuzzyASAExample.c ExampleCost.c asa.c ------------8<------------ bottom cut <- top ------------->8------------ ============================================================================== ======================================================================== @@CHIPSET for ASA_QUEUE 3 Nov 96 chipset is public domain source code by Paul Sander for doubly-linked lists, from ftp.cc.utexas.edu:/source/utilities/misc/chipset/. This is an alternative source of implementing lists for ASA_QUEUE to that given below by Helmut Jarausch. ======================================================================== @@Hashed Lists for ASA_QUEUE 15 Dec 95 Helmut Jarausch wrote: Dear Lester, I have tried to downdate my previous C++ version to C. Since I have no access to a pre-ANSI C-compiler I could test the code only with an ANSI C compiler. It accepts the code with HAVE_ANSI either TRUE or FALSE. The results are the same. If you like to include this code into ASA, there can be a further simplification by fixing (and if you like, inlining) the functions fct = cost_function Raster -> inline Hash -> probably inline Probably you have noticed that the approximation with my code is (mathematically) slightly different from your code. To make the difference clear it suffices to look at a single real variable since in multi-D each dimension is treated the same. Both approaches approximate a function by a piecewise constant function. Let Queue_Resolution = 0.5. Given the following random sequence out of the interval [0,10] 3.96 8.40 3.53 4.45 3.19 8.86 0.16 5.84 1.59 3.84 My code returns fct at these pts 4.0 8.5 3.5 4.5 3.0 9.0 0.0 6.0 1.5 4.0 The current ASA returns fct at these pts 3.96 8.40 3.96 3.96 3.19 8.40 0.16 5.84 1.59 3.96 I think both schemes are reasonable. For my scheme I need an "offset" value. The most reasonable suggestions are either Minimum or (Minimum+Maximum)/2. The formula which can be used in the function Raster is round( (x-x_off) / resolution ) * resolution Of course a user could define a fancy resolution scheme doing a nonlinear gridding like taking logarithms, etc. But I do need the value of x_off. A fixed prescribed value will not do. Assume Minimum = 1E50, Maximum = 1E51, and assume x_off has been set to 0. Then there is danger that round((x-x_off)/resolution) would result in an integer bigger than INT_MAX, in which case the result is set to INT_MAX. Thank you for your attention, best wishes, Helmut. /* --8<--------8<---- CUT HERE TO CUT-LINE BELOW ----8<--------8<-- */ #include #include #include #ifdef HAVE_ANSI #define Const const #else #define Const #endif #define Local static typedef double RET_TYPE; typedef double ARG_TYPE; typedef RET_TYPE (*FCT) (Const ARG_TYPE *); typedef void (*RASTER_FCT) (int, ARG_TYPE *, Const ARG_TYPE *); typedef unsigned (*HASH_FCT) (int, Const ARG_TYPE *); struct FCE { RET_TYPE fct_val; ARG_TYPE *arg; struct FCE *Hash_next, *prev, *next; }; Local struct { int Dim_Arg; unsigned Size; int Eval_Count, Eval_Hit, Saved_Arg_Index; struct FCE **Hash_Table, *Queue, *Head, *Tail; ARG_TYPE *Saved_Arg, *Temp_Arg; FCT fct; RASTER_FCT Raster; HASH_FCT Hash; } A_F_C; /* The constructor initializes an object. Here it take the size of the cache, the dimension of the user function's argument, the name of the user function itself, a function (Raster) (see below) and the primitive Hash function. */ #ifdef HAVE_ANSI Local int A_F_C_Ctor (int Cache_Size, int Arg_dim, FCT fct_fct, RASTER_FCT Raster_fct, HASH_FCT Hash_fct) #else Local int A_F_C_Ctor (Cache_Size, Arg_dim, fct_fct, Raster_fct, Hash_fct) int Cache_Size, Arg_dim; FCT fct_fct; RASTER_FCT Raster_fct; HASH_FCT Hash_fct; #endif { struct FCE *Q; int i; A_F_C.Dim_Arg = Arg_dim; A_F_C.fct = fct_fct; A_F_C.Raster = Raster_fct; A_F_C.Hash = Hash_fct; /* compute smallest prime number not less than Cache_Size */ { int q, p = Cache_Size; if (!(p & 1)) p++; p -= 2; do { p += 2; q = 3; while (p % q && q * q <= p) q += 2; } while (q * q <= p); A_F_C.Size = p; } A_F_C.Eval_Count = A_F_C.Eval_Hit = 0; A_F_C.Saved_Arg_Index = A_F_C.Size; A_F_C.Hash_Table = (struct FCE **) malloc (sizeof (struct FCE *) * A_F_C.Size); if (!A_F_C.Hash_Table) return FALSE; A_F_C.Queue = (struct FCE *) malloc (sizeof (struct FCE) * A_F_C.Size); if (!A_F_C.Queue) return FALSE; A_F_C.Saved_Arg = (ARG_TYPE *) malloc (sizeof (ARG_TYPE) * A_F_C.Dim_Arg * A_F_C.Size); if (!A_F_C.Saved_Arg) return FALSE; A_F_C.Temp_Arg = (ARG_TYPE *) malloc (sizeof (ARG_TYPE) * A_F_C.Dim_Arg); if (!A_F_C.Temp_Arg) return FALSE; Q = A_F_C.Head = A_F_C.Queue; for (i = 0; i < A_F_C.Size; i++) A_F_C.Hash_Table[i] = 0; for (i = 0; i < A_F_C.Size - 1; i++) { Q->next = Q + 1; (Q + 1)->prev = Q; Q++; } A_F_C.Tail = Q; A_F_C.Tail->next = A_F_C.Head; A_F_C.Head->prev = A_F_C.Tail; return TRUE; } Local void A_F_C_Dtor () { free (A_F_C.Temp_Arg); free (A_F_C.Saved_Arg); free (A_F_C.Queue); free (A_F_C.Hash_Table); } /* updating the MRU queue */ #ifdef HAVE_ANSI Local void A_F_C_requeue (struct FCE *Q) #else Local void A_F_C_requeue (Q) struct FCE *Q; #endif { /* dequeue Q */ Q->prev->next = Q->next; Q->next->prev = Q->prev; if (Q == A_F_C.Tail) A_F_C.Tail = Q->prev; /* make Q the new Head */ Q->next = A_F_C.Head; Q->prev = A_F_C.Head->prev; Q->prev->next = Q; A_F_C.Head->prev = Q; A_F_C.Head = Q; } /* locate an entry in the hash list */ #ifdef HAVE_ANSI Local void A_F_C_Hash_find (Const ARG_TYPE * Parm, struct FCE **QP, unsigned *HV) #else Local void A_F_C_Hash_find (Parm, QP, HV) ARG_TYPE *Parm; struct FCE **QP; unsigned *HV; #endif { int i; struct FCE *Q; *HV = A_F_C.Hash (A_F_C.Dim_Arg, Parm) % A_F_C.Size; Q = A_F_C.Hash_Table[*HV]; while (Q) /* typically this loop will take only 1 or 2 times */ { i = 0; while (i < A_F_C.Dim_Arg && Parm[i] == (Q->arg)[i]) i++; if (i == A_F_C.Dim_Arg) break; /* found */ Q = Q->Hash_next; } *QP = Q; } /* to remove an entry that is to be phased out */ #ifdef HAVE_ANSI Local void A_F_C_Hash_delete (struct FCE *Q_old) #else Local void A_F_C_Hash_delete (Q_old) struct FCE *Q_old; #endif { int i; unsigned HV; struct FCE *Father, *Q; HV = A_F_C.Hash (A_F_C.Dim_Arg, Q_old->arg) % A_F_C.Size; Q = A_F_C.Hash_Table[HV]; Father = 0; while (Q) { i = 0; while (i < A_F_C.Dim_Arg && (Q_old->arg)[i] == (Q->arg)[i]) i++; if (i == A_F_C.Dim_Arg) break; /* found */ Father = Q; Q = Q->Hash_next; } if (!Q) exit (13); /* Hash delete internal error - should never occur */ if (Father) Father->Hash_next = Q->Hash_next; else A_F_C.Hash_Table[HV] = Q->Hash_next; } #ifdef HAVE_ANSI Local RET_TYPE A_F_C_Eval (Const ARG_TYPE * Parm) #else Local RET_TYPE A_F_C_Eval (Parm) ARG_TYPE *Parm; #endif { int i; unsigned HV; struct FCE *Q; A_F_C.Raster (A_F_C.Dim_Arg, A_F_C.Temp_Arg, Parm); A_F_C.Eval_Count++; A_F_C_Hash_find (A_F_C.Temp_Arg, &Q, &HV); if (Q) /* found */ { A_F_C.Eval_Hit++; A_F_C_requeue (Q); /* move to the Head */ } else /* not found */ { Q = A_F_C.Tail; /* use least recently used entry */ A_F_C_requeue (Q); if (A_F_C.Saved_Arg_Index) /* queue not yet filled */ Q->arg = &A_F_C.Saved_Arg[--A_F_C.Saved_Arg_Index * A_F_C.Dim_Arg]; else /* delete old value from A_F_C.Hash_Table */ A_F_C_Hash_delete (Q); for (i = 0; i < A_F_C.Dim_Arg; i++) (Q->arg)[i] = A_F_C.Temp_Arg[i]; Q->fct_val = A_F_C.fct (A_F_C.Temp_Arg); /* compute function fresh */ /* insert a new entry into the hash list */ Q->Hash_next = A_F_C.Hash_Table[HV]; A_F_C.Hash_Table[HV] = Q; } return Q->fct_val; } Local int A_F_C_Hit_Ratio () { return ((100 * A_F_C.Eval_Hit) / A_F_C.Eval_Count); } /* ---------------- Just an example for checking it ------------------------- */ double Res[] = {0.2, 1.0, 2.0}; /* this corresponds to ASA's Queue_Resolution */ #define round(x) ( x < 0 ? -(int)(0.5-x) : (int)(x+0.5) ) #ifdef HAVE_ANSI void raster (int N, double Y[], Const double X[]) #else void raster (N, Y, X) int N; double *X, *Y; #endif { int i; for (i = 0; i < N; i++) Y[i] = round (X[i] / Res[i]) * Res[i]; } #ifdef HAVE_ANSI unsigned Hash (int N, Const double X[]) #else unsigned Hash (N, X) int N; double *X; #endif { unsigned H = 0; int i; for (i = 0; i < N; i++) H ^= round (X[i] / Res[i]); return H; } #ifdef HAVE_ANSI double myfct (Const double X[]) #else double myfct (X) double *X; #endif { double F = 0; int i; for (i = 0; i < 3; i++) F += fabs (X[i]) / Res[i]; return F; } int main () { int Cache_Sz, i, j, N; double X[3]; double Buffer[200][3]; double diff, diff_max, FE, mF; Cache_Sz = 200; /* create and initialize */ if (!A_F_C_Ctor (Cache_Sz, 3, myfct, raster, Hash)) { exit (1); } printf ("Res[1:3]: "); scanf ("%lg %lg %lg", Res, Res + 1, Res + 2); printf ("current resolutions are %g %g %g\n", Res[0], Res[1], Res[2]); diff_max = 0; for (i = 0; i < 200; i++) { for (j = 0; j < 3; j++) X[j] = Buffer[i][j] = drand48 () * 10.0; FE = A_F_C_Eval (X); mF = myfct (X); diff = fabs (FE - mF); if (diff > diff_max) diff_max = diff; } printf ("Hit_ratio: %d\n", A_F_C_Hit_Ratio ()); /* Now replay all the old values - this should result a hit ratio of 100 per cent for this part */ for (i = 0; i < 200; i++) { for (j = 0; j < 3; j++) X[j] = Buffer[i][j]; A_F_C_Eval (X); } printf ("Hit_ratio: %d\n", A_F_C_Hit_Ratio ()); /* free storage */ A_F_C_Dtor (); return 0; } /* -->8-------->8---- CUT HERE TO CUT-LINE ABOVE ---->8-------->8-- */ ======================================================================== @@FDLIBM Code for FDLIBM_EXP, FDLIBM_LOG and FDLIBM_POW This is an example of FDLIBM code described in the ASA-README, to replace the standard exp(), log(), and/or pow() functions, which may yield some speed-up in ASA performance. This code should first be tested with the standard ASA_TEST OPTIONS to see if the resulting asa_out file agrees with the test_asa file. For a copy of FDLIBM, send a message "send index from fdlibm" to netlib@research.att.com, or retrieve via FTP from netlib.att.com:/netlib/fdlibm The code below can be appended to to user.c (as of version 12.3) or to asa.c (as of version 12.2). FDLIBM_EXP, FDLIBM_LOG and/or FDLIBM_POW can be defined as TRUE in your complation procedures. /* --8<--------8<---- CUT HERE TO CUT-LINE BELOW ----8<--------8<-- */ #if FDLIBM_POW #define FDLIBM TRUE #else #if FDLIBM_LOG #define FDLIBM TRUE #else #if FDLIBM_EXP #define FDLIBM TRUE #endif #endif #endif #ifndef FDLIBM #define FDLIBM FALSE #endif #if FDLIBM #ifdef __LITTLE_ENDIAN #define __HI(x) *(1+(int*)&x) #define __LO(x) *(int*)&x #define __HIp(x) *(1+(int*)x) #define __LOp(x) *(int*)x #else #define __HI(x) *(int*)&x #define __LO(x) *(1+(int*)&x) #define __HIp(x) *(int*)x #define __LOp(x) *(1+(int*)x) #endif #endif /* These defines block out some original code. */ #ifndef FDLIBM_POW_ORG #define FDLIBM_POW_ORG FALSE #endif #ifndef FDLIBM_LOG_ORG #define FDLIBM_LOG_ORG FALSE #endif #ifndef FDLIBM_EXP_ORG #define FDLIBM_EXP_ORG FALSE #endif #if FDLIBM #if HAVE_ANSI static const double #else static double #endif huge = 1.0e+300; #endif /* FDLIBM */ #if FDLIBM_POW /* @(#)e_pow.c 1.2 95/01/04 */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== */ /* __ieee754_pow(x,y) return x**y * n * Method: Let x = 2 * (1+f) * 1. Compute and return log2(x) in two pieces: * log2(x) = w1 + w2, * where w1 has 53-24 = 29 bit trailing zeros. * 2. Perform y*log2(x) = n+y' by simulating muti-precision * arithmetic, where |y'|<=0.5. * 3. Return x**y = 2**n*exp(y'*log2) * * Special cases: * 1. (anything) ** 0 is 1 * 2. (anything) ** 1 is itself * 3. (anything) ** NAN is NAN * 4. NAN ** (anything except 0) is NAN * 5. +-(|x| > 1) ** +INF is +INF * 6. +-(|x| > 1) ** -INF is +0 * 7. +-(|x| < 1) ** +INF is +0 * 8. +-(|x| < 1) ** -INF is +INF * 9. +-1 ** +-INF is NAN * 10. +0 ** (+anything except 0, NAN) is +0 * 11. -0 ** (+anything except 0, NAN, odd integer) is +0 * 12. +0 ** (-anything except 0, NAN) is +INF * 13. -0 ** (-anything except 0, NAN, odd integer) is +INF * 14. -0 ** (odd integer) = -( +0 ** (odd integer) ) * 15. +INF ** (+anything except 0,NAN) is +INF * 16. +INF ** (-anything except 0,NAN) is +0 * 17. -INF ** (anything) = -0 ** (-anything) * 18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer) * 19. (-anything except 0 and inf) ** (non-integer) is NAN * * Accuracy: * pow(x,y) returns x**y nearly rounded. In particular * pow(integer,integer) * always returns the correct integer provided it is * representable. * * Constants : * The hexadecimal values are the intended ones for the following * constants. The decimal values may be used, provided that the * compiler will convert from decimal to binary accurately enough * to produce the hexadecimal values shown. */ #if FDLIBM_POW_ORG #include "fdlibm.h" #ifdef __STDC__ static const double #else static double #endif #else /* FDLIBM_POW_ORG */ #if HAVE_ANSI static const double #else static double #endif #endif /* FDLIBM_POW_ORG */ bp[] = {1.0, 1.5,}, dp_h[] = {0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */ dp_l[] = {0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */ #if FDLIBM_POW_ORG zero = 0.0, one = 1.0, two = 2.0, two53 = 9007199254740992.0, /* 0x43400000, 0x00000000 */ huge = 1.0e300, tiny = 1.0e-300, /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */ #endif /* FDLIBM_POW_ORG */ two53 = 9007199254740992.0, /* 0x43400000, 0x00000000 */ tiny = 1.0e-300, /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */ L1 = 5.99999999999994648725e-01, /* 0x3FE33333, 0x33333303 */ L2 = 4.28571428578550184252e-01, /* 0x3FDB6DB6, 0xDB6FABFF */ L3 = 3.33333329818377432918e-01, /* 0x3FD55555, 0x518F264D */ L4 = 2.72728123808534006489e-01, /* 0x3FD17460, 0xA91D4101 */ L5 = 2.30660745775561754067e-01, /* 0x3FCD864A, 0x93C9DB65 */ L6 = 2.06975017800338417784e-01, /* 0x3FCA7E28, 0x4A454EEF */ P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */ P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */ P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */ P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */ P5 = 4.13813679705723846039e-08, /* 0x3E663769, 0x72BEA4D0 */ lg2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */ lg2_h = 6.93147182464599609375e-01, /* 0x3FE62E43, 0x00000000 */ lg2_l = -1.90465429995776804525e-09, /* 0xBE205C61, 0x0CA86C39 */ ovt = 8.0085662595372944372e-0017, /* -(1024-log2(ovfl+.5ulp)) */ cp = 9.61796693925975554329e-01, /* 0x3FEEC709, 0xDC3A03FD =2/(3ln2) */ cp_h = 9.61796700954437255859e-01, /* 0x3FEEC709, 0xE0000000 =(float)cp */ cp_l = -7.02846165095275826516e-09, /* 0xBE3E2FE0, 0x145B01F5 =tail of cp_h */ ivln2 = 1.44269504088896338700e+00, /* 0x3FF71547, 0x652B82FE =1/ln2 */ ivln2_h = 1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2 */ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail */ #if FDLIBM_POW_ORG #ifdef __STDC__ double __ieee754_pow (double x, double y) #else double __ieee754_pow (x, y) double x, y; #endif #else /* FDLIBM_POW_ORG */ #if HAVE_ANSI double s_pow (double x, double y) #else double s_pow (x, y) double x, y; #endif #endif /* FDLIBM_POW_ORG */ { double z, ax, z_h, z_l, p_h, p_l; double y1, t1, t2, r, s, t, u, v, w; int i, j, k, yisint, n; int hx, hy, ix, iy; unsigned lx, ly; hx = __HI (x); lx = __LO (x); hy = __HI (y); ly = __LO (y); ix = hx & 0x7fffffff; iy = hy & 0x7fffffff; /* y==zero: x**0 = 1 */ if ((iy | ly) == 0) return 1.0; /* +-NaN return x+y */ if (ix > 0x7ff00000 || ((ix == 0x7ff00000) && (lx != 0)) || iy > 0x7ff00000 || ((iy == 0x7ff00000) && (ly != 0))) return x + y; /* determine if y is an odd int when x < 0 * yisint = 0 ... y is not an integer * yisint = 1 ... y is an odd int * yisint = 2 ... y is an even int */ yisint = 0; if (hx < 0) { if (iy >= 0x43400000) yisint = 2; /* even integer y */ else if (iy >= 0x3ff00000) { k = (iy >> 20) - 0x3ff; /* exponent */ if (k > 20) { j = ly >> (52 - k); if ((j << (52 - k)) == ly) yisint = 2 - (j & 1); } else if (ly == 0) { j = iy >> (20 - k); if ((j << (20 - k)) == iy) yisint = 2 - (j & 1); } } } /* special value of y */ if (ly == 0) { if (iy == 0x7ff00000) { /* y is +-inf */ if (((ix - 0x3ff00000) | lx) == 0) return y - y; /* inf**+-1 is NaN */ else if (ix >= 0x3ff00000) /* (|x|>1)**+-inf = inf,0 */ return (hy >= 0) ? y : 0.0; else /* (|x|<1)**-,+inf = inf,0 */ return (hy < 0) ? -y : 0.0; } if (iy == 0x3ff00000) { /* y is +-1 */ if (hy < 0) return 1.0 / x; else return x; } if (hy == 0x40000000) return x * x; /* y is 2 */ if (hy == 0x3fe00000) { /* y is 0.5 */ if (hx >= 0) /* x >= +0 */ return sqrt (x); } } ax = fabs (x); /* special value of x */ if (lx == 0) { if (ix == 0x7ff00000 || ix == 0 || ix == 0x3ff00000) { z = ax; /*x is +-0,+-inf,+-1 */ if (hy < 0) z = 1.0 / z; /* z = (1/|x|) */ if (hx < 0) { if (((ix - 0x3ff00000) | yisint) == 0) { z = (z - z) / (z - z); /* (-1)**non-int is NaN */ } else if (yisint == 1) z = -z; /* (x<0)**odd = -(|x|**odd) */ } return z; } } /* (x<0)**(non-int) is NaN */ if ((((hx >> 31) + 1) | yisint) == 0) return (x - x) / (x - x); /* |y| is huge */ if (iy > 0x41e00000) { /* if |y| > 2**31 */ if (iy > 0x43f00000) { /* if |y| > 2**64, must o/uflow */ if (ix <= 0x3fefffff) return (hy < 0) ? huge * huge : tiny * tiny; if (ix >= 0x3ff00000) return (hy > 0) ? huge * huge : tiny * tiny; } /* over/underflow if x is not close to one */ if (ix < 0x3fefffff) return (hy < 0) ? huge * huge : tiny * tiny; if (ix > 0x3ff00000) return (hy > 0) ? huge * huge : tiny * tiny; /* now |1-x| is tiny <= 2**-20, suffice to compute log(x) by x-x^2/2+x^3/3-x^4/4 */ t = x - 1; /* t has 20 trailing zeros */ w = (t * t) * (0.5 - t * (0.3333333333333333333333 - t * 0.25)); u = ivln2_h * t; /* ivln2_h has 21 sig. bits */ v = t * ivln2_l - w * ivln2; t1 = u + v; __LO (t1) = 0; t2 = v - (t1 - u); } else { double s2, s_h, s_l, t_h, t_l; n = 0; /* take care subnormal number */ if (ix < 0x00100000) { ax *= two53; n -= 53; ix = __HI (ax); } n += ((ix) >> 20) - 0x3ff; j = ix & 0x000fffff; /* determine interval */ ix = j | 0x3ff00000; /* normalize ix */ if (j <= 0x3988E) k = 0; /* |x|> 1) | 0x20000000) + 0x00080000 + (k << 18); t_l = ax - (t_h - bp[k]); s_l = v * ((u - s_h * t_h) - s_h * t_l); /* compute log(ax) */ s2 = s * s; r = s2 * s2 * (L1 + s2 * (L2 + s2 * (L3 + s2 * (L4 + s2 * (L5 + s2 * L6))))); r += s_l * (s_h + s); s2 = s_h * s_h; t_h = 3.0 + s2 + r; __LO (t_h) = 0; t_l = r - ((t_h - 3.0) - s2); /* u+v = s*(1+...) */ u = s_h * t_h; v = s_l * t_h + t_l * s; /* 2/(3log2)*(s+...) */ p_h = u + v; __LO (p_h) = 0; p_l = v - (p_h - u); z_h = cp_h * p_h; /* cp_h+cp_l = 2/(3*log2) */ z_l = cp_l * p_h + p_l * cp + dp_l[k]; /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */ t = (double) n; t1 = (((z_h + z_l) + dp_h[k]) + t); __LO (t1) = 0; t2 = z_l - (((t1 - t) - dp_h[k]) - z_h); } s = 1.0; /* s (sign of result -ve**odd) = -1 else = 1 */ if ((((hx >> 31) + 1) | (yisint - 1)) == 0) s = -1.0; /* (-ve)**(odd int) */ /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */ y1 = y; __LO (y1) = 0; p_l = (y - y1) * t1 + y * t2; p_h = y1 * t1; z = p_l + p_h; j = __HI (z); i = __LO (z); if (j >= 0x40900000) { /* z >= 1024 */ if (((j - 0x40900000) | i) != 0) /* if z > 1024 */ return s * huge * huge; /* overflow */ else { if (p_l + ovt > z - p_h) return s * huge * huge; /* overflow */ } } else if ((j & 0x7fffffff) >= 0x4090cc00) { /* z <= -1075 */ if (((j - 0xc090cc00) | i) != 0) /* z < -1075 */ return s * tiny * tiny; /* underflow */ else { if (p_l <= z - p_h) return s * tiny * tiny; /* underflow */ } } /* * compute 2**(p_h+p_l) */ i = j & 0x7fffffff; k = (i >> 20) - 0x3ff; n = 0; if (i > 0x3fe00000) { /* if |z| > 0.5, set n = [z+0.5] */ n = j + (0x00100000 >> (k + 1)); k = ((n & 0x7fffffff) >> 20) - 0x3ff; /* new k for n */ t = 0.0; __HI (t) = (n & ~(0x000fffff >> k)); n = ((n & 0x000fffff) | 0x00100000) >> (20 - k); if (j < 0) n = -n; p_h -= t; } t = p_l + p_h; __LO (t) = 0; u = t * lg2_h; v = (p_l - (t - p_h)) * lg2 + t * lg2_l; z = u + v; w = v - (z - u); t = z * z; t1 = z - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5)))); r = (z * t1) / (t1 - 2.0) - (w + z * w); z = 1.0 - (r - z); j = __HI (z); j += (n << 20); #if FDLIBM_POW_ORG if ((j >> 20) <= 0) z = scalbn (z, n); /* subnormal output */ else __HI (z) += (n << 20); #else /* FDLIBM_POW_ORG */ __HI (z) += (n << 20); #endif /* FDLIBM_POW_ORG */ return s * z; } #endif /* FDLIBM_POW */ #if FDLIBM_LOG /* @(#)e_log.c 1.2 95/01/04 */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== */ /* __ieee754_log(x) * Return the logarithm of x * * Method : * 1. Argument Reduction: find k and f such that * x = 2^k * (1+f), * where sqrt(2)/2 < 1+f < sqrt(2) . * * 2. Approximation of log(1+f). * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) * = 2s + 2/3 s**3 + 2/5 s**5 + ....., * = 2s + s*R * We use a special Reme algorithm on [0,0.1716] to generate * a polynomial of degree 14 to approximate R The maximum error * of this polynomial approximation is bounded by 2**-58.45. In * other words, * 2 4 6 8 10 12 14 * R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s * (the values of Lg1 to Lg7 are listed in the program) * and * | 2 14 | -58.45 * | Lg1*s +...+Lg7*s - R(z) | <= 2 * | | * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2. * In order to guarantee error in log below 1ulp, we compute log * by * log(1+f) = f - s*(f - R) (if f is not too large) * log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy) * * 3. Finally, log(x) = k*ln2 + log(1+f). * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo))) * Here ln2 is split into two floating point number: * ln2_hi + ln2_lo, * where n*ln2_hi is always exact for |n| < 2000. * * Special cases: * log(x) is NaN with signal if x < 0 (including -INF) ; * log(+INF) is +INF; log(0) is -INF with signal; * log(NaN) is that NaN with no signal. * * Accuracy: * according to an error analysis, the error is always less than * 1 ulp (unit in the last place). * * Constants: * The hexadecimal values are the intended ones for the following * constants. The decimal values may be used, provided that the * compiler will convert from decimal to binary accurately enough * to produce the hexadecimal values shown. */ #if FDLIBM_LOG_ORG #include "fdlibm.h" #ifdef __STDC__ static const double #else static double #endif #else /* FDLIBM_LOG_ORG */ #if HAVE_ANSI static const double #else static double #endif #endif /* FDLIBM_LOG_ORG */ ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */ ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */ two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */ Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */ Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */ Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */ Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */ Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */ Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */ Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */ #if FDLIBM_LOG_ORG static double zero = 0.0; #ifdef __STDC__ double __ieee754_log (double x) #else double __ieee754_log (x) double x; #endif #else /* FDLIBM_LOG_ORG */ #if HAVE_ANSI double s_log (double x) #else double s_log (x) double x; #endif #endif /* FDLIBM_LOG_ORG */ { double hfsq, f, s, z, R, w, t1, t2, dk; int k, hx, i, j; unsigned lx; hx = __HI (x); /* high word of x */ lx = __LO (x); /* low word of x */ k = 0; if (hx < 0x00100000) { /* x < 2**-1022 */ if (((hx & 0x7fffffff) | lx) == 0) #if FDLIBM_LOG_ORG return -two54 / zero; /* log(+-0)=-inf */ #else return -huge; /* log(+-0)=-inf */ #endif if (hx < 0) return (x - x) / 0.0; /* log(-#) = NaN */ k -= 54; x *= two54; /* subnormal number, scale up x */ hx = __HI (x); /* high word of x */ } if (hx >= 0x7ff00000) return x + x; k += (hx >> 20) - 1023; hx &= 0x000fffff; i = (hx + 0x95f64) & 0x100000; __HI (x) = hx | (i ^ 0x3ff00000); /* normalize x or x/2 */ k += (i >> 20); f = x - 1.0; if ((0x000fffff & (2 + hx)) < 3) { /* |f| < 2**-20 */ if (f == 0.0) if (k == 0) return 0.0; else { dk = (double) k; return dk * ln2_hi + dk * ln2_lo; } R = f * f * (0.5 - 0.33333333333333333 * f); if (k == 0) return f - R; else { dk = (double) k; return dk * ln2_hi - ((R - dk * ln2_lo) - f); } } s = f / (2.0 + f); dk = (double) k; z = s * s; i = hx - 0x6147a; w = z * z; j = 0x6b851 - hx; t1 = w * (Lg2 + w * (Lg4 + w * Lg6)); t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7))); i |= j; R = t2 + t1; if (i > 0) { hfsq = 0.5 * f * f; if (k == 0) return f - (hfsq - s * (hfsq + R)); else return dk * ln2_hi - ((hfsq - (s * (hfsq + R) + dk * ln2_lo)) - f); } else { if (k == 0) return f - s * (f - R); else return dk * ln2_hi - ((s * (f - R) - dk * ln2_lo) - f); } } #endif /* FDLIBM_LOG */ #if FDLIBM_EXP /* @(#)e_exp.c 1.2 95/01/04 */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== */ /* __ieee754_exp(x) * Returns the exponential of x. * * Method * 1. Argument reduction: * Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658. * Given x, find r and integer k such that * * x = k*ln2 + r, |r| <= 0.5*ln2. * * Here r will be represented as r = hi-lo for better * accuracy. * * 2. Approximation of exp(r) by a special rational function on * the interval [0,0.34658]: * Write * R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ... * We use a special Reme algorithm on [0,0.34658] to generate * a polynomial of degree 5 to approximate R. The maximum error * of this polynomial approximation is bounded by 2**-59. In * other words, * R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5 * (where z=r*r, and the values of P1 to P5 are listed below) * and * | 5 | -59 * | 2.0+P1*z+...+P5*z - R(z) | <= 2 * | | * The computation of exp(r) thus becomes * 2*r * exp(r) = 1 + ------- * R - r * r*R1(r) * = 1 + r + ----------- (for better accuracy) * 2 - R1(r) * where * 2 4 10 * R1(r) = r - (P1*r + P2*r + ... + P5*r ). * * 3. Scale back to obtain exp(x): * From step 1, we have * exp(x) = 2^k * exp(r) * * Special cases: * exp(INF) is INF, exp(NaN) is NaN; * exp(-INF) is 0, and * for finite argument, only exp(0)=1 is exact. * * Accuracy: * according to an error analysis, the error is always less than * 1 ulp (unit in the last place). * * Misc. info. * For IEEE double * if x > 7.09782712893383973096e+02 then exp(x) overflow * if x < -7.45133219101941108420e+02 then exp(x) underflow * * Constants: * The hexadecimal values are the intended ones for the following * constants. The decimal values may be used, provided that the * compiler will convert from decimal to binary accurately enough * to produce the hexadecimal values shown. */ #if FDLIBM_EXP_ORG #include "fdlibm.h" #ifdef __STDC__ static const double #else static double #endif #else /* FDLIBM_EXP_ORG */ #if HAVE_ANSI static const double #else static double #endif #endif /* FDLIBM_EXP_ORG */ halF[2] = {0.5, -0.5,}, #if FDLIBM_EXP_ORG one = 1.0, huge = 1.0e+300, #endif twom1000 = 9.33263618503218878990e-302, /* 2**-1000=0x01700000,0 */ o_threshold = 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */ u_threshold = -7.45133219101941108420e+02, /* 0xc0874910, 0xD52D3051 */ ln2HI[2] = {6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */ -6.93147180369123816490e-01,}, /* 0xbfe62e42, 0xfee00000 */ ln2LO[2] = {1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */ -1.90821492927058770002e-10,}, /* 0xbdea39ef, 0x35793c76 */ #if FDLIBM_POW invln2 = 1.44269504088896338700e+00; /* 0x3ff71547, 0x652b82fe */ #else invln2 = 1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */ P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */ P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */ P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */ P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */ P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */ #endif #if FDLIBM_EXP_ORG #ifdef __STDC__ double __ieee754_exp (double x) /* default IEEE double exp */ #else double __ieee754_exp (x) /* default IEEE double exp */ double x; #endif #else /* FDLIBM_EXP_ORG */ #if HAVE_ANSI double s_exp (double x) /* default IEEE double exp */ #else double s_exp (x) /* default IEEE double exp */ double x; #endif #endif /* FDLIBM_EXP_ORG */ { double y, hi, lo, c, t; int k, xsb; unsigned hx; #if FDLIBM_EXP_ORG #else /* turn off warnings that these may be used uninitialized */ hi = lo = 0.0; k = 0; #endif hx = __HI (x); /* high word of x */ xsb = (hx >> 31) & 1; /* sign bit of x */ hx &= 0x7fffffff; /* high word of |x| */ /* filter out non-finite argument */ if (hx >= 0x40862E42) { /* if |x|>=709.78... */ if (hx >= 0x7ff00000) { if (((hx & 0xfffff) | __LO (x)) != 0) return x + x; /* NaN */ else return (xsb == 0) ? x : 0.0; /* exp(+-inf)={inf,0} */ } if (x > o_threshold) return huge * huge; /* overflow */ if (x < u_threshold) return twom1000 * twom1000; /* underflow */ } /* argument reduction */ if (hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */ if (hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */ hi = x - ln2HI[xsb]; lo = ln2LO[xsb]; k = 1 - xsb - xsb; } else { k = invln2 * x + halF[xsb]; t = k; hi = x - t * ln2HI[0]; /* t*ln2HI is exact here */ lo = t * ln2LO[0]; } x = hi - lo; } else if (hx < 0x3e300000) { /* when |x|<2**-28 */ if (huge + x > 1.0) return 1.0 + x; /* trigger inexact */ } else k = 0; /* x is now in primary range */ t = x * x; c = x - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5)))); if (k == 0) return 1.0 - ((x * c) / (c - 2.0) - x); else y = 1.0 - ((lo - (x * c) / (2.0 - c)) - hi); if (k >= -1021) { __HI (y) += (k << 20); /* add k to y's exponent */ return y; } else { __HI (y) += ((k + 1000) << 20); /* add k to y's exponent */ return y * twom1000; } } #endif /* FDLIBM_EXP */ /* -->8-------->8---- CUT HERE TO CUT-LINE ABOVE ---->8-------->8-- */ ======================================================================== ======================================================================== @@RATFOR vfsr.r Code (1987) unsupported Note that this example of VFSR was set up to fit a Lagrangian corresponding to a set of coupled stochastic differential equations. Any other function cost() could be used instead, and the main program vfsr and the common block vfsr_com.r need only set up any data requirements for that problem. Then, it likely is necessary to include the subroutines trajec and stat, and any of the erf...() functions (now a standard call in most math libraries). Also note that I wrote ASA for a reason: ASA is better than VFSR in all ways, except possibly ease of use if you cannot get by using the easy entry into ASA described in the ASA-README file. ----------8<---------- top cut 1/2 -> bottom 1/2 ----------->8---------- # RATFOR vfsr.r for JANUS project (1987) # Requires common block vfsr_com.r (1987) # # ===== program vfsr { include "vfsr_com.r"; # IMSL ctime() also in prtint open (unit=uout,file='runvfsr_dat',status='new'); call cputim; open (unit=uin,file='vfsrdat_dat',status='old'); # CASES: program vfsr, function cost, file vfsrdat_dat write (uout,*) 'Lanchester drift, additive noise'; # write (uout,*) 'Lanchester drift, multiplicative noise'; # write (uout,*) 'Pade drift, additive noise'; # write (uout,*) 'Pade drift, multiplicative noise'; # write (uout,*) 'Pade-Lanchester drift, additive noise'; onbc = 0; write (uout,*) ' '; read (uin,*) (dname(ndum),ndum=1,2); write (uout,*) (dname(ndum),ndum=1,2); read (uin,*) pvar,ktraj; write (uout,*) 'pvar=',pvar,'ktraj=',ktraj; read (uin,*) (dname(ndum),ndum=1,2); write (uout,*) (dname(ndum),ndum=1,2); read (uin,*) tbeg,tfin; write (uout,*) 'tbeg=',tbeg,'tfin=',tfin; write (uout,*) ' '; read (uin,*) (dname(ndum),ndum=1,2); write (uout,*) (dname(ndum),ndum=1,2); read (uin,*) cycle,seed; write (uout,*) 'cycle=',cycle,'seed=',seed; write (uout,*) ' '; read (uin,*) (dname(ndum),ndum=1,3); write (uout,*) (dname(ndum),ndum=1,3); read (uin,*) instop,rstop,cfac; # rstop alternative to instop write (uout,*) 'instop=',instop,'rstop=',rstop,'cfac=',cfac; write (uout,*) ' '; pi1 = dacos(- one); pi2 = two*pi1; expfl = five*five; logln = dlog(ten); pardim = 0; write (uout,*) ' '; read (uin,*) (dname(ndum),ndum=1,5); write (uout,*) (dname(ndum),ndum=1,5); do p = 1,pvar { read (uin,*) ndum,mpar(p),sys(p),varmin(p),varmax(p); mlag(p) = mpar(p) - sys(p); write (uout,*) 'var_p=',ndum,'mpar(p)=',mpar(p),'sys(p)=',sys(p),'mlag(p)=',mlag(p),'varmin=',varmin(p),'varmax=',varmax(p); pardim = pardim + mpar(p); } write (uout,*) ' '; read (uin,*) (dname(ndum),ndum=1,3); write (uout,*) (dname(ndum),ndum=1,3); do p = 1,pvar { do m = 1,mpar(p) { read (uin,*) par(p,m),parmin(p,m),parmax(p,m); begpar(p,m) = par(p,m); pargen(p,m) = par(p,m); # set tempc0's tempc0(p,m) = one; write (uout,*) '(p,m)=',p,m,':',par(p,m),parmin(p,m),parmax(p,m); } } write (uout,*) ' '; close (uin); # warm up uniform random number generator do ndum = 1,1000 { dum1 = unif (); } stop; call trajec; close (utraj); # statistics on trajectories action = cost(); deltat = dfloat(ktraj)*(tfin - tbeg)/(dfloat(npnt1)); begact = action; finact = action; finlag = lagst; findet = logdet; # set tempb0 tempb0 = dabs(action); tmpbsv = tempb0; # nstop alternative to instop dum1 = one; nstop = int(action/dum1); nstop = instop; write (uout,*) 'tempb0=initial_cost=',action,'nstop=',nstop,'rstop=',rstop; write (uout,*) 'initial_lagst=',lagst,'initial_logdet=',logdet; write (uout,*) ' '; write (uout,*) '====='; do ndum = 1,cycle { tempb0 = tmpbsv; write (uout,*) ' seed=',seed; call anneal; write (uout,*) ' '; write (uout,*) '====='; write (uout,*) 'results_for_cycle=',ndum; write (uout,*) ' rstop=',rstop,'ratio=',ratio; write (uout,*) ' instop=',instop,'nstop=',nstop,'nacc=',nacc; dum1 = dfloat(nacc)/dfloat(ngen); write (uout,*) 'nacc/ngen=',dum1,'tempb=',tempb,'cost=',action; do p = 1,pvar { do m = 1,mpar(p) { write (uout,*) '(p,m)=',p,m,'par=',par(p,m),'tempc=',tempc(p,m); } } write (uout,*) '_____'; write (uout,*) 'min_cost=',endact,'tot_gen=',totgen,'endtb=',endtb; action = endact; do p = 1,pvar { do m = 1,mpar(p) { pargen(p,m) = endpar(p,m); } } call senstv; do p = 1,pvar { do m = 1,mpar(p) { write (uout,*) '(p,m)=',p,m,'min_par=',endpar(p,m),'sens=',sens(p,m),'curv=',curv(p,m); } } if (endact < finact) { finact = endact; finlag = endlag; findet = enddet; do p = 1,pvar { do m = 1,mpar(p) { finpar(p,m) = endpar(p,m); } } } write (uout,*) '====='; } call stat; write (uout,*) 'final_cost=',finact,'final_lagst=',finlag,'deltat=',deltat,'findet=',findet; do p = 1,pvar { do m = 1,mpar(p) { write (uout,*) '(p,m)=',p,m,'finpar=',finpar(p,m); } } close (uout); stop; end; } # ===== # ===== subroutine trajec { include "vfsr_com.r"; npoint = 0; npnt1 = 0; open (unit=utraj,file='vfsrtrj_dat',status='old'); tmax = 0; write (uout,*) ' '; read (utraj,*) (dname(ndum),ndum=1,1); write (uout,*) (dname(ndum),ndum=1,1); read (utraj,*) (dname(ndum),ndum=1,1),((ddname(p,y),y=1,sys(p)),p=1,pvar); write (uout,*) (dname(ndum),ndum=1,1),((ddname(p,y),y=1,sys(p)),p=1,pvar); do k = 1,ktraj { read (utraj,*) sktime(k); write (uout,*) ' trajk=',k,'sktime_init=',sktime(k); sktim2(k) = 0; s2 = 0; do s = 1,sktime(k) { read (utraj,*) dum1,((ddum(p,y),y=1,sys(p)),p=1,pvar); if (s == 1) { do p = 1,pvar { do y = 1,sys(p) { var0(p,y) = ddum(p,y); } } write (uout,*) 'var0:',((var0(p,y),y=1,sys(p)),p=1,pvar); } if (dum1 >= tbeg & dum1 <= tfin) { sktim2(k) = sktim2(k) + 1; s2 = s2 + 1; time(s2,k) = dum1; do p = 1,pvar { do y = 1,sys(p) { var(p,y,s2,k) = ddum(p,y); # use fractional forces # var(p,y,s2,k) = ddum(p,y)/var0(p,y); # } } } } sktime(k) = sktim2(k); write (uout,*) 'trajk=',k,'sktime_qual=',sktime(k); if (sktime(k) > tmax) { tmax= sktime(k); } npoint = npoint + sktime(k); npnt1 = npnt1 + sktime(k) - 1; } write (uout,*) ' '; write (uout,*) 'npoint=',npoint; return; end; } # ===== # ===== subroutine stat { include "vfsr_com.r"; write (uout,*) '=-=-='; write (uout,*) ' ','time','av1','var1','av2','var2'; # check increments of time tim = time(1,1) - (time(2,1) - time(1,1)); # check increments of time do s = 1,tmax { tim = tim + (time(2,1) - time(1,1)); dum1 = 0; dum2 = 0; dum3 = 0; dum4 = 0; ndum = 0; y = 1; do k = 1,ktraj { if (s <= sktime(k)){ ndum = ndum + 1; dum1 = dum1 + var(1,y,s,k); dum2 = dum2 + var(2,y,s,k); dum3 = dum3 + var(1,y,s,k)**2; dum4 = dum4 + var(2,y,s,k)**2; } } dum1 = dum1/ndum; dum2 = dum2/ndum; dum3 = dum3/ndum; dum4 = dum4/ndum; dum3 = dsqrt(dum3 - dum1**2); dum4 = dsqrt(dum4 - dum2**2); rdum1 = real(dum1); rdum2 = real(dum2); rdum3 = real(dum3); rdum4 = real(dum4); write (uout,*) tim,rdum1,rdum3,rdum2,rdum4; } write (uout,*) '=-=-='; write (uout,*) 'time','mean(1)','var(1)','mean(2)','var(2)'; ndum = sktime(1) - 1; do s = 1,ndum { do p = 1,pvar { q(p) = zero; dq(p) = zero; } do k = 1,ktraj { do p = 1,pvar { do y = 1,sys(p) { yp = y + mlag(p); q(p) = q(p) + pargen(p,yp)*var(p,y,s,k); dq(p) = dq(p) + (pargen(p,yp)*var(p,y,s,k))**2; } } } do p = 1,pvar { q(p) = q(p)/ktraj; dq(p) = dq(p)/ktraj - q(p)**2; } write (uout,*) time(s,1),(q(p),dq(p),p=1,pvar); } return; end; } # ===== # ===== subroutine anneal { include "vfsr_com.r"; totgen = 0; intgen = 0; intacc = 0; ratio = dfloat(intacc + 1)/dfloat(intgen + 1); ngen = 1; nacc = 1; naccsv = nacc; naccb = 1; tempb = tempb0; endtb = tempb0; savact = begact; savlag = lagst; savdet = logdet; endact = action; endlag = lagst; enddet = logdet; # temperature Tf=T0*10**-dum1 at time tf=10**dum2 dum1 = two*two; dum2 = two*two; dum3 = dum1*logln; dum4 = dum2*logln; cexp = dum3*dexp(-dum4/dfloat(pardim)); write (uout,*) 'cexp=',cexp; do p = 1,pvar { do m = 1,mpar(p) { endpar(p,m) = begpar(p,m); vfrgen(p,m) = 1; } } while ( (nacc <= nstop) & (nacc < instop) & (ratio > rstop) ) { dum1 = dexp(-logln*two*two*three); do p = 1,pvar { do m = 1,mpar(p) { dum3 = -cexp*dfloat(vfrgen(p,m))**(one/dfloat(pardim)); dum2 = expchk(dum3); tempc(p,m) = tempc0(p,m)*dexp(dum2); if (tempc(p,m) < dum1) { write (uout,*) '*** tempc too small','(p,m)=',p,m,tempc(p,m); return; } } } dum1 = dexp(-logln*ten); dum3 = -cexp*dfloat(naccb)**(one/dfloat(pardim)); dum2 = expchk(dum3); tempb = tempb0*dexp(dum2); if (tempb < dum1) { write (uout,*) '*** tempb too small',tempb; return; } call gener; action = cost(); call accept; ratio = dfloat(intacc + 1)/dfloat(intgen + 1); if (action < endact) { intacc = 0; intgen = 0; endact = action; endlag = lagst; enddet = logdet; endtb = tempb; do p = 1,pvar { do m = 1,mpar(p) { endpar(p,m) = pargen(p,m); } } } ndum = mod(nacc,cfac); if ( (ndum == 0) & (naccsv == nacc) ) { call senstv; call prtint; # re-anneal dum3 = zero; dum4 = one; if (dum4 > dum3) { do p = 1,pvar { do m = 1,mpar(p) { if (parmin(p,m) != parmax(p,m)) { dum1 = tempc(p,m)*(snsmax/sens(p,m)); if (dum1 < tempc0(p,m)) { vfrgen(p,m) = int(((dlog(tempc0(p,m)/dum1))/cexp)**pardim); } else { vfrgen(p,m) = 1; } } } } if (tempb0 > dabs(savact)) { tempb0 = dabs(savact); } if (tempb > dabs(endact)) { naccb = int(((dlog(tempb0/dabs(endact)))/cexp)**pardim); } else { naccb = int(((dlog(tempb0/tempb))/cexp)**pardim); } } } } return; end; } # ===== # ===== subroutine prtint { include "vfsr_com.r"; write (uout,*) '_____'; write (uout,*) 'intermediate_results:'; call cputim; write (uout,*) 'savact=',savact,'ngen=',ngen,'nacc=',nacc,'ratio=',ratio; write (uout,*) 'endact=',endact,'naccb=',naccb,'tempb0=',tempb0,'tempb=',tempb; write (uout,*) ' savlag=',savlag,'endlag=',endlag,'savdet=',savdet,'enddet=',enddet,'deltat=',deltat; write (uout,*) '(p,m):','endpar=','tempc=','vfrgen=','curv=','sens='; do p = 1,pvar { do m = 1,mpar(p) { write (uout,*) p,m,':',endpar(p,m),tempc(p,m),vfrgen(p,m),curv(p,m),sens(p,m); } } return; end; } # ===== # ===== function cost() { include "vfsr_com.r"; cost = zero; lagst = zero; logdet = zero; do k = 1,ktraj { ndum = sktime(k) - 1; do s = 1,ndum { dt = time(s+1,k) - time(s,k); do p = 1,pvar { q(p) = zero; qp(p) = zero; } do p = 1,pvar { do y = 1,sys(p) { yp = y + mlag(p); q(p) = q(p) + pargen(p,yp)*var(p,y,s,k); qp(p) = qp(p) + pargen(p,yp)*var(p,y,s+1,k); } dq(p) = qp(p) - q(p); } # DRIFTS # Lanchester drift: drift(1) = pargen(1,1)*q(2) + pargen(1,2)*q(1)*q(2); drift(2) = pargen(2,1)*q(1) + pargen(2,2)*q(2)*q(1); # # Lanchester time-dependent drift: # drift(1) = pargen(1,1)*q(2) + pargen(1,2)*q(1)*q(2) + pargen(1,4)*q(2)*time(s,k); # drift(2) = pargen(2,1)*q(1) + pargen(2,2)*q(2)*q(1) + pargen(1,4)*q(1)*time(s,k); # # Pade drift: # drift(1) = (pargen(1,1)*q(2) + pargen(1,2)*q(1)*q(2)) / ( one + pargen(1,4)*q(1) + pargen(1,5)*q(2) ); # drift(2) = (pargen(2,1)*q(1) + pargen(2,2)*q(2)*q(1)) / ( one + pargen(2,4)*q(2) + pargen(2,5)*q(1) ); # # Pade-Lanchester drift: # drift(1) = (pargen(1,4)*q(2) + pargen(1,5)*q(1)*q(2)) / ( one + pargen(1,1)*q(1) + pargen(1,2)*q(2) ); # drift(2) = (pargen(2,4)*q(1) + pargen(2,5)*q(2)*q(1)) / ( one + pargen(2,1)*q(2) + pargen(2,2)*q(1) ); # # DIFFUSIONS # multiplicative noise: # diffus(1,1) = ( pargen(1,3)*(one+q(1)) )**2; # diffus(2,2) = ( pargen(2,3)*(one+q(2)) )**2; # # additive noise: diffus(1,1) = pargen(1,3)**2; diffus(2,2) = pargen(2,3)**2; # diffus(1,2) = zero; diffus(2,1) = diffus(1,2); dum1 = one/(diffus(1,1)*diffus(2,2)); if (diffus(1,2) != zero) { dum1 = dum1 - one/diffus(1,2)**2; } detq = one/dum1; do p = 1,pvar { hq(p) = dq(p) - drift(p)*dt; } do p = 1,pvar { lagp = hq(p)**2/(two*diffus(p,p)*dt); lagst = lagst + drift(p)**2/(two*diffus(p,p)); if (onbc > 0) { # add b.c. if (qp(p) - varmin(p) < two * dsqrt(diffus(p,p)*dt)) { dum1 = hq(p) + two*(qp(p) - varmin(p)); dum2 = (dum1**2 + four*drift(p)*dt*(q(p) - varmin(p)))/(two*diffus(p,p)*dt); dum3 = expchk(-dum2); lagp = - dlog(dexp(-lagp) - dexp(-dum3)); } else if (varmax(p) - qp(p) < two * dsqrt(diffus(p,p)*dt)) { dum1 = qp(p) - two*varmax(p) + q(p) + drift(p)*dt; dum2 = dum1 - two*drift(p)*dt; dum1 = erf(dum1/dsqrt(two*diffus(p,p)*dt)); dum2 = (dum2**2 + four*drift(p)*dt*(q(p) - varmax(p)))/(two*diffus(p,p)*dt); dum2 = expchk(-dum2); dum3 = drift(p)/diffus(p,p); dum4 = expchk(two*dum3*(qp(p) - varmax(p))); dum1 = dum1*dexp(expchk(two*dum3*(qp(p) - varmax(p)))); dum3 = dum3*dsqrt(pi2*dt*diffus(p,p)); dum1 = dum3*(dum1 + dexp(dum4)) + dexp(-dum2); lagp = - dlog(dexp(-lagp) + dum1); } # end b.c. } cost = cost + lagp; do pp = 1,pvar { if (pp > p & diffus(p,pp) != zero) { cost = cost + hq(p)*hq(pp)/(diffus(p,pp)*dt); lagst = lagst + drift(p)*drift(pp)/diffus(p,pp); } } } logdet = dfloat(pvar)*dlog(pi2*dt)/two + dlog(detq)/two; cost = cost + logdet; } } cost = cost / dfloat(npnt1); lagst = lagst / dfloat(npnt1); logdet = logdet / dfloat(npnt1); return; end; } # ===== # ===== subroutine gener { include "vfsr_com.r"; do p = 1,pvar { ppvar = p; do m = 1,mpar(p) { mmpar = m; repeat { dum1 = par(p,m) + vryfst()*(parmax(p,m) - parmin(p,m)); totgen = totgen + 1; } until (dum1 <= parmax(p,m) & dum1 >= parmin(p,m)); pargen(p,m) = dum1; vfrgen(p,m) = vfrgen(p,m) + 1; } } intgen = intgen + 1; ngen = ngen + 1; return; end; } # ===== # ===== subroutine accept { include "vfsr_com.r"; dact = (action - savact)/tempb; naccsv = naccsv + 1; if (dact > expfl) { dum1 = zero; dum2 = one; } else if (dact < - expfl) { dum1 = one; dum2 = zero; } else { dum1 = dexp(-dact); dum2 = unif (); } if (dum1 > dum2) { savact = action; savlag = lagst; savdet = logdet; do p = 1,pvar { do m = 1,mpar(p) { par(p,m) = pargen(p,m); } } intacc = intacc + 1; nacc = nacc + 1; naccsv = nacc; naccb = naccb + 1; } return; end; } # ===== # ===== function vryfst() { include "vfsr_com.r"; # runs at about 1000 vryfst's/sec on VAX-785 # significant digit accuracy = log10(# generated points) dum1 = unif (); dum2 = tempc(ppvar,mmpar); if (dum1 < half) { dum3 = -one; } else { dum3 = one; } vryfst = dum3*dum2*((one + one/dum2)**dabs(two*dum1 - one) - one); return; end; } # ===== # ===== function unif () { include "vfsr_com.r"; data d2p31 /2147483648.0d0/,d2p31m /2147483647.0d0/,dmult /16807.0d0/; seed = dmod(dmult*seed,d2p31m); unif = seed/d2p31; return; end; } # ===== # ===== subroutine senstv { include "vfsr_com.r"; do p = 1,pvar { do m = 1,mpar(p) { pargen(p,m) = endpar(p,m); } } do p = 1,pvar { do m = 1,mpar(p) { if (parmin(p,m) != parmax(p,m)) { dpar = 0.01d0*endpar(p,m); pargen(p,m) = 1.01d0*endpar(p,m); dum1 = cost(); pargen(p,m) = 0.99d0*endpar(p,m); dum2 = cost(); sens(p,m) = dabs((dum1 - endact)/dpar)*(parmax(p,m) - parmin(p,m)); curv(p,m) = (dum1 - two*endact + dum2)/(endpar(p,m)*dpar**2); pargen(p,m) = endpar(p,m); } } } snsmax = sens(1,1); do p = 1,pvar { do m = 1,mpar(p) { if (parmin(p,m) != parmax(p,m)) { if (sens(p,m) > snsmax) { snsmax = sens(p,m); } } } } return; end; } # ===== # ===== subroutine cputim { include "vfsr_com.r"; # IMSL ctime() # dum1 = ctime(); # write (uout,*) 'cputime(sec)=',dum1; return; end; } # ===== # ===== function expchk(dum1) { include "vfsr_com.r"; if (dum1 > expfl) { expchk = expfl; } else if (dum1 < -expfl) { expchk = -expfl; } else { expchk = dum1; } return; end; } # ===== # ===== function erf(x) { integer sn; double precision erf, erf4, erfc4, x; if (x < 0.0) { sn = -1; x = -x; } else sn = 1; if (x > 8.0) erf = 1.0; else if (x < 1.5) erf = erf4(x); else erf = 1.0 - erfc4(x); erf = erf * sn; return; end; } # ===== # ===== function erf4(x) # erf good for small x { double precision erf4, x, x2, sum, sqrtpi; double precision t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12; data t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12/_ 0.6666667,0.2666667,0.07619048,0.01693122,3.078403e-3,_ 4.736005e-4,6.314673e-5,7.429027e-6,7.820028e-7,_ 7.447646e-8,6.476214e-9/; x2 = x*x; sum = t5+x2*(t6+x2*(t7+x2*(t8+x2*(t9+x2*(t10+x2*(t11+x2*t12)))))); sqrtpi = dsqrt(dacos(-1.0d0)); erf4 = 2.0*exp(-x2)/sqrtpi*(x*(1+x2*(t2+x2*(t3+x2*(t4+x2*sum))))); return; end; } # ===== # ===== function erfc4(x) # erfc good for large x { double precision erfc4, x, x2, sum, v, sqrtpi; x2 = x*x; v = 1.0/(2.0*x2); sum = v/(1+8*v/(1+9*v/(1+10*v/(1+11*v/(1+12*v))))); sum = v/(1+3*v/(1+4*v/(1+5*v/(1+6*v/(1+7*sum))))); sqrtpi = dsqrt(dacos(-1.0d0)); erfc4 = 1.0/(exp(x2)*x*sqrtpi*(1+v/(1+2*sum))); return; end; } # ===== # ===== function erfc(x) { double precision erfc, erf4, erfc4, erf, x; if (x < 0.0) { erfc = 1.0 - erf(x); return; } if ( x > 8.0) erfc = 0.0; else if (x < 1.5) erfc = 1.0 - erf4(x); else erfc = erfc4(x); return; end; } # ===== ----------8<---------- bottom cut 1/2 <- top 1/2 ----------->8---------- RATFOR vfsr_com.r Common Block (1987) unsupported ----------8<---------- top cut 2/2 -> bottom 2/2 ----------->8---------- # RATFOR vfsr_com.r for JANUS project (1987) # Common block for vfsr.r (1987) # #----- # vfsr_com.r character*8 dname(10),ddname(10,10); real rdum1,rdum2,rdum3,rdum4; integer k,m,p,pp,s,s2,ndum,mpar,pvar,mlag,sys,y,yp,mmpar,ppvar,ktraj,cycle,npoint,npnt1,sktime,sktim2,tmax,onbc,_ instop,nstop,ngen,nacc,naccsv,naccb,intgen,intacc,uin,udata,utraj,uout,totgen,cfac,vfrgen,pardim; double precision dum1,dum2,dum3,dum4,ddum(2,10),cost,expchk,tim,dt,deltat,detq,pi1,pi2,expfl,action,begact,savact,dact,endact,finact,endtb,_ varmin,varmax,lagst,savlag,endlag,lagp,finlag,logdet,savdet,enddet,findet,_ tempb,tempb0,tmpbsv,par,parmin,parmax,pargen,begpar,endpar,finpar,_ var,var0,time,q(2),qp(2),dq(2),hq(2),pref(2),erf,_ drift,diffus,tempc,tempc0,cauchy,vryfst,unif,gaussr,seed,d2p31,d2p31m,dmult,_ ratio,rstop,sens,curv,snsmax,dpar,zero,one,two,three,four,five,seven,ten,half,cexp,logln,tbeg,tfin; parameter (uin=14,utraj=15,uout=16,udata=17); common /sa/ pvar,mpar(2),mlag(2),sys(2),mmpar,ppvar,ktraj,cycle,npoint,npnt1,sktime(100),sktim2(50),tim,deltat,tmax,seed,instop,nstop,totgen,pardim,_ cfac,vfrgen(2,10),ngen,nacc,naccsv,naccb,intgen,intacc,pi1,pi2,expfl,_ action,begact,savact,endact,finact,endtb,tempb,tbeg,tfin,_ drift(2),diffus(2,2),varmin(5),varmax(5),lagst,savlag,endlag,finlag,logdet,savdet,enddet,findet,_ tempb0,par(2,10),parmin(2,10),parmax(2,10),pargen(2,10),begpar(2,10),endpar(2,10),finpar(2,10),_ var(2,10,200,50),var0(2,10),time(200,50),tempc(2,10),tempc0(2,10),ratio,rstop,sens(2,10),curv(2,10),snsmax,cexp,logln,onbc; external ctime; data zero/0.0d0/,one/1.0d0/,two/2.0d0/,three/3.0d0/,four/4.0d0/,five/5.0d0/,seven/7.0d0/,ten/10.0d0/,half/0.5d0/; #----- ----------8<---------- bottom cut 2/2 <- top 2/2 ----------->8---------- ======================================================================== Lester Ingber Copyright (c) 1994-2008 Lester Ingber. All Rights Reserved, Unless Stated Otherwise. $Id: asa_contrib.txt,v 10.46 2008/02/19 00:32:26 ingber Exp ingber $