https://www.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 https://www.ingber.com/asa_examples.txt https://www.ingber.com/asa_examples.txt For papers dealing with difficult systems across a wide range of disciplines, see https://www.ingber.com/asa_papers.html https://www.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 @@Python and OpenOpt 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. These can be accessed from https://www.ingber.com/private/ASA_PARALLEL_BINDER_2001_01_23.tar.gz 01 Jan 14 When ASA_PARALLEL is set to TRUE, code in ASA-MAKEFILE, asa.c, asa_usr.c, and asa_usr_asa.h is set to run under OpenMP using the gcc compiler; search for instances of _OPENMP. Change/add lines of code to correspond to your choice of parallel algorithm and compiler. ======================================================================== @@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 16 Mar 12 [updated] http://ssakata.sdf.org/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 https://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 . ============================================================================== @@Python and OpenOpt Interface with ASA 14 Jul 12 Robert Jordens has written a Python interface to ASA at http://pypi.python.org/pypi/pyasa A OpenOpt solver using this Python interface is being developed by Dmitrey Kroshko at http://openopt.org/GLP ============================================================================== @@FUZZY_ASA 5 Mar 11 As of ASA-28.2, a ASA_FUZZY code is part of ASA; it has been modified to make it easier for users to select ASA OPTIONS used for tuning. For more examples, and current ASA_FUZZY code, go the the developer's URL below. Caution: Often ASA_FUZZY will turn on QUENCHing > 1, violating the proof of ASA. For many systems, this speeding up of the sampling process can be a welcome efficiency, but in some systems global minima may be missed. 5 May 09 Hime Junior has contributed asa_fuzzy.c to use with ASA. 20 Jun 10 Hime has created a site with his asa_fuzzy code and also a Windows graphical interface: http://sites.google.com/site/stochasticglobaloptimization/home/fuzzy-asa The code follows these guidelines: In Makefile: Add asa_fuzzy.o to USEROBJS In Makefile or asa_usr_asa.h: Add -DASA_FUZZY=TRUE to DEFINE_OPTIONS or put "#define ASA_FUZZY TRUE" into asa_usr_asa.h Similarly, be sure QUENCH_PARAMETERS and QUENCH_COST are set to TRUE asa_usr.h Add: #if ASA_FUZZY int InitFuzzyASA (int NoOfDimensions); void CloseFuzzyASA (void); void FuzzyControl (double *x, double fvalue, USER_DEFINES * USER_OPTIONS, int NoParam); #endif /* ASA_FUZZY */ asa_usr.c: Be sure initial values of QUENCH_PARAMETERS and QUENCH_COST are set at their proper sites. Just before any call to asa (...) put: #if ASA_FUZZY InitFuzzyASA (*parameter_dimension); #endif /* ASA_FUZZY */ Just after any call to asa() put: #if ASA_FUZZY CloseFuzzyASA (); #endif /* ASA_FUZZY */ In cost_function(...): Just after the value of the cost is calculated (say value is costXX) and *cost_flag is set, put: #if ASA_FUZZY if (*cost_flag == TRUE && USER_OPTIONS->Locate_Cost == 12) FuzzyControl (x, costXX, USER_OPTIONS, *parameter_dimension); #endif /* ASA_FUZZY */ Using Locate_Cost can be useful to skip transient phases, calls during numerical derivatives, etc. In the above example, Locate_Cost is used to permit use of FuzzyControl() when "the cost function is being used for a new generated state just after a new best state was achieved" (see ASA-README.[]). In asa_fuzzy.c you may decide to change the number of samples of the cost-function values used for the fuzzy "energy" evaluation, set at #define NOOFSAMPLES 30 ------------8<------------ top cut -> bottom ------------->8------------ /* asa_fuzzy.c contributed by Hime Junior 20090505 */ #include "asa_usr.h" #define THRESHOLD1 -0.0001 #define THRESHOLDDEVIATION 5e-4 #define NOOFSAMPLES 30 #define PERFORMANCE_TARGET 0.1 #define FACTOR_a 2 /* As NOOFSAMPLES regulates the size of arrays storing data about some best values found recently, the user can freely change it. Please, notice that values higher than 100 or lower than 5, for instance, could cause problems to the fuzzy controller's "reasoning". */ static double *FuzzyParameters[NOOFSAMPLES]; static double FuzzyValues[NOOFSAMPLES + 1]; static double FuzzyMinima[NOOFSAMPLES + 1]; /* Initialization of fuzzy controller's data structures - should be called just before activations of asa() */ extern int InitFuzzyASA (int NoOfDimensions) { int index; for (index = 0; index < NOOFSAMPLES; ++index) if ((FuzzyParameters[index] = (double *) calloc (NoOfDimensions, sizeof (double))) == NULL) return -1; return 0; } /* Release of fuzzy controller's data structures obtained previously - should be called just after returning of asa() */ extern void CloseFuzzyASA (void) { int index; for (index = 0; index < NOOFSAMPLES; ++index) free (FuzzyParameters[index]); } void MeanAndDeviation (double *Vector, int NumberOfElements, double *Mean, double *Deviation) { int i; double aux; *Mean = *Deviation = 0; for (i = 0; i < NumberOfElements; i++) { aux = Vector[i]; *Mean += aux; *Deviation += aux * aux; } *Mean /= NumberOfElements; *Deviation /= NumberOfElements; *Deviation = *Deviation - (*Mean) * (*Mean); *Deviation = sqrt (fabs (*Deviation)); } double SubEnergy (double InputValue, double Minimum) { double valuhat, argulog; valuhat = InputValue - Minimum; argulog = 1 / (1 + exp (-(valuhat + FACTOR_a))); if (argulog == 0) argulog = 1e-100; return log (argulog); } double DeltaFactor (double MeanSub) { #define DELTATPOS .3 #define DELTATNULO 0 #define DELTATNEG -.3 #define VALMINLOC log(1/(1+exp(-FACTOR_a))) #define INTER1 .1 double MembershipMeanZero = 0; double MembershipMeanMedium = 0; double ResuCrisp; /* RULE #1 - IF MeanSub is ZERO THEN decreasing rate is POSITIVE ( We are distant from present basic local minimum ) */ if (MeanSub > THRESHOLD1) MembershipMeanZero = (MeanSub - THRESHOLD1) / (-THRESHOLD1); /* RULE #2 - IF MeanSub is VALMINLOC ( FUZZY NUMBER ) THEN decreasing rate is POSITIVE ( We are in a region near the present basic local minimum ) */ if (MeanSub >= VALMINLOC && MeanSub <= VALMINLOC + INTER1) MembershipMeanMedium = (VALMINLOC + INTER1 - MeanSub) / INTER1; ResuCrisp = MembershipMeanZero * DELTATPOS + MembershipMeanMedium * DELTATPOS; return ResuCrisp; } void AlterQuench (USER_DEFINES * PresentASAOptions, int NoParam, double Mean, double Deviation) { int i, j; double Delta, Meanaux, Deviationaux; double aux[NOOFSAMPLES]; Delta = DeltaFactor (Mean); if (PresentASAOptions->User_Quench_Cost_Scale[0] < 100) PresentASAOptions->User_Quench_Cost_Scale[0] *= (1 + .5 * Delta); for (i = 0; i < NoParam; i++) if (PresentASAOptions->User_Quench_Param_Scale[i] < 100) PresentASAOptions->User_Quench_Param_Scale[i] *= (1 + .2 * Delta); if (Deviation < THRESHOLDDEVIATION) { for (i = 0; i < NoParam; i++) if (PresentASAOptions->User_Quench_Param_Scale[i] < 100 && PresentASAOptions->User_Quench_Param_Scale[i] > .1) { for (j = 0; j < NOOFSAMPLES; j++) aux[j] = FuzzyParameters[j][i]; MeanAndDeviation (aux, NOOFSAMPLES, &Meanaux, &Deviationaux); PresentASAOptions->User_Quench_Param_Scale[i] /= (1 + .15 * Delta * exp (-Deviationaux)); } } } extern void FuzzyControl (double *x, double fvalue, USER_DEFINES * USER_OPTIONS, int NoParam) { static double ActualPerformance, Mean, Deviation, aux[NOOFSAMPLES]; static int IndVal = 0; int i; IndVal++; if (IndVal % (NOOFSAMPLES + 1)) { FuzzyValues[IndVal] = fvalue; FuzzyMinima[IndVal] = *USER_OPTIONS->Best_Cost; for (i = 0; i < NoParam; i++) FuzzyParameters[IndVal - 1][i] = USER_OPTIONS->Best_Parameters[i]; /* Stores better results until now */ return; } IndVal = 0; if (FuzzyMinima[1] != 0) ActualPerformance = (FuzzyMinima[1] - FuzzyMinima[NOOFSAMPLES]) / fabs (FuzzyMinima[1]); else return; if (ActualPerformance > PERFORMANCE_TARGET) return; for (i = 0; i < NOOFSAMPLES; i++) aux[i] = SubEnergy (FuzzyValues[i + 1], FuzzyMinima[i + 1]); // Zero based MeanAndDeviation (aux, NOOFSAMPLES, &Mean, &Deviation); AlterQuench (USER_OPTIONS, NoParam, Mean, Deviation); return; } ------------8<------------ bottom cut <- top ------------->8------------ 2 May 09 I have added to asa_paper.html %A Hime Aguiar e Oliveira Jr. %T Fuzzy Modeling with Adaptive Simulated Annealing %C Rio de Janeiro, Brazil %D 2004 %O URL http://www.optimization-online.org/DB_HTML/2004/07/901.html and two other 2007 and 2009 papers. 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 above 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. ============================================================================== ======================================================================== @@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-2016 Lester Ingber. All Rights Reserved, Unless Stated Otherwise. $Id: asa_contrib.txt,v 10.69 2016/12/25 22:53:51 ingber Exp ingber $