ExternalVendorCode/trick
1
0
Fork 0
mirror of https://github.com/nasa/trick.git synced 2024-12-25 07:41:08 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Delete unused baseball directory in trick models.

Browse Source
This commit is contained in:
John M. Penn 2015-08-31 15:01:44 -05:00
parent c98a662ae6
commit 77c9293c29
21 changed files with 0 additions and 1320 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
trick_models/baseball/aero/include/cannon_aero.d
View File

@ -1,24 +0,0 @@
/* Initialize cannon ball shot */
CANNON_AERO.lift_method = Smits_Smith ;
CANNON_AERO.coefficient_drag = 0.45 ;
CANNON_AERO.coefficient_cross = 0.044 ;
CANNON_AERO.mass = 0.145 ;
CANNON_AERO.air_density = 1.29 ;
CANNON_AERO.ball_radius {cm} = 3.63 ;
CANNON_AERO.ball_area {cm2} = 41.59 ;
CANNON_AERO.g = -9.81 ;
CANNON_AERO.force_jet_Z_plus {lbf} = 2.0 ;
/* Regula Falsi impact critter setup */
#define BIG_TGO 10000
CANNON_AERO.rf.lower_set = No ;
CANNON_AERO.rf.upper_set = No ;
CANNON_AERO.rf.iterations = 0 ;
CANNON_AERO.rf.fires = 0 ;
CANNON_AERO.rf.x_lower = BIG_TGO ;
CANNON_AERO.rf.t_lower = BIG_TGO ;
CANNON_AERO.rf.x_upper = BIG_TGO ;
CANNON_AERO.rf.t_upper = BIG_TGO ;
CANNON_AERO.rf.delta_time = BIG_TGO ;
CANNON_AERO.rf.error_tol = 1.0e-9 ;
CANNON_AERO.rf.mode = Decreasing ;

71
trick_models/baseball/aero/include/cannon_aero.h
View File

@ -1,71 +0,0 @@
/*********************************** TRICK HEADER **************************
PURPOSE: (Test Baseball)
***************************************************************************/
#ifndef CANNON_AERO_H
#define CANNON_AERO_H
#include "sim_services/include/regula_falsi.h"
typedef enum {
Hard_Coded_Coefficient_Lift, /* You come up with Cl */
Smits_Smith, /* Cl ~= 0.54*S^0.4 (S = rw/V) */
Adair_Giordano, /* Lift_Force = mass*0.00041*w_cross_V */
Tombras /* Cl = 1/(2.022 + 0.981v/w) */
} Lift_Estimation_Method ;
typedef struct {
double pos[3] ; /* M position */
double vel[3] ; /* M/s velocity */
double acc[3] ; /* M/s2 acceleration */
double omega[3] ; /* r/s Angular velocity of cannonball */
double theta ; /* r Angle from x-axis to axis rotation */
double phi ; /* r Angle from z-axis to axis rotation */
double omega0 ; /* r/s Magnitude of angular velocity about
axis of rotation */
/* Impact */
REGULA_FALSI rf ; /* -- Dynamic event params for impact */
int impact ; /* -- Has impact occured */
double impact_pos ; /* M How far ball lands in field */
/* Forces */
double force_gravity[3] ; /* N Gravitational force */
double force_drag[3] ; /* N Drag force opposite dir of velocity */
double force_magnus[3] ; /* N Force due to spin, dir grav cross drag */
double force_cross[3] ; /* N Side directional force */
double force_total[3] ; /* N Sum of all forces */
/* Force magnitudes */
double mag_force_drag ; /* N Magnitude of drag force */
double mag_force_magnus ; /* N Magnitude of magnus force */
double mag_force_cross ; /* N Magnitude of cross force */
double mag_omega ; /* r/s Magnitude of angular velocity */
void** force_collect ; /* -- For collect statement */
/* Environment and Properties */
double mass ; /* kg Mass of cannonball */
double air_density ; /* kg/M3 Air density at 20C */
double ball_radius ; /* M Radius of cannonball */
double ball_area ; /* M2 Cross sectional area of ball */
double spin_parameter ; /* -- S=r*omega/speed */
double g ; /* M/s2 Gravitational acceleration */
/* Coefficients drag, lift and cross */
Lift_Estimation_Method lift_method ; /* -- How to find lift force */
double coefficient_drag ; /* -- Drag coefficient */
double coefficient_lift ; /* -- Lift coefficient */
double coefficient_cross ; /* -- Cross-Force coefficient */
/* Jet */
int jet_on ; /* -- 0|1 */
int jet_count ; /* -- How many jet firings? */
double force_jet[3] ; /* N Jet force per firing */
double force_jet_Z_plus ; /* N Configurable force of jet in Z+ direction */
double time_to_fire_jet_1 ; /* s First jet firing time */
double time_to_fire_jet_2 ; /* s Second jet firing time */
double time_to_fire_jet_3 ; /* s Third jet firing time */
double time_to_fire_jet_4 ; /* s Fourth jet firing time */
} CANNON_AERO ;
#endif

24
trick_models/baseball/aero/include/cannon_integ_aero.d
View File

@ -1,24 +0,0 @@
#define NUM_STEP 12
#define NUM_VARIABLES 6
INTEGRATOR.state = alloc(NUM_VARIABLES) ;
INTEGRATOR.deriv = alloc(NUM_STEP) ;
INTEGRATOR.state_ws = alloc(NUM_STEP) ;
for (int kk = 0 ; kk < NUM_STEP ; kk++ ) {
INTEGRATOR.deriv[kk] = alloc(NUM_VARIABLES) ;
INTEGRATOR.state_ws[kk] = alloc(NUM_VARIABLES) ;
}
INTEGRATOR.stored_data = alloc(8) ;
for (int kk = 0 ; kk < 8 ; kk++ ) {
INTEGRATOR.stored_data[kk] = alloc(NUM_VARIABLES) ;
}
INTEGRATOR.option = Runge_Kutta_4 ;
INTEGRATOR.init = True ;
INTEGRATOR.first_step_deriv = Yes ;
INTEGRATOR.num_state = NUM_VARIABLES ;
#undef NUM_STEP
#undef NUM_VARIABLES

31
trick_models/baseball/aero/src/cannon_collect_forces.c
View File

@ -1,31 +0,0 @@
/*********************************** TRICK HEADER **************************
PURPOSE: (Collect all forces and calculate acceleration)
***************************************************************************/
#include "../include/cannon_aero.h"
#include "sim_services/include/collect_macros.h"
int cannon_collect_forces(
CANNON_AERO *C )
{
double **collected_forces ;
int ii ;
/* Collect external forces on the ball */
collected_forces = (double**)(C->force_collect) ;
C->force_total[0] = 0.0 ;
C->force_total[1] = 0.0 ;
C->force_total[2] = 0.0 ;
for( ii = 0 ; ii < NUM_COLLECT(collected_forces) ; ii++ ) {
C->force_total[0] += collected_forces[ii][0] ;
C->force_total[1] += collected_forces[ii][1] ;
C->force_total[2] += collected_forces[ii][2] ;
}
/* Solve for xyz acceleration */
C->acc[0] = C->force_total[0] / C->mass ;
C->acc[1] = C->force_total[1] / C->mass ;
C->acc[2] = C->force_total[2] / C->mass ;
return 0 ;
}

27
trick_models/baseball/aero/src/cannon_force_cross.c
View File

@ -1,27 +0,0 @@
/*********************************** TRICK HEADER **************************
PURPOSE: (Cross Force or Side Force )
***************************************************************************/
#include "../include/cannon_aero.h"
#include "trick_utils/math/include/trick_math.h"
int cannon_force_cross(
CANNON_AERO *C )
{
double magnus_cross_drag[3] ;
double norm_magnus_cross_drag[3] ;
double k, speed ;
/* k = 1/2*rho*Cy*A*V^2 */
speed = V_MAG( C->vel ) ;
k = (-0.5)*C->air_density*C->coefficient_cross*C->ball_area*speed*speed ;
/* F = k*(M x D)/|M x D| */
V_CROSS( magnus_cross_drag, C->force_magnus, C->force_drag ) ;
V_NORM( norm_magnus_cross_drag, magnus_cross_drag ) ;
V_SCALE( C->force_cross, norm_magnus_cross_drag, k) ;
C->mag_force_cross = V_MAG( C->force_cross ) ;
return 0 ;
}

24
trick_models/baseball/aero/src/cannon_force_drag.c
View File

@ -1,24 +0,0 @@
/*********************************** TRICK HEADER **************************
PURPOSE: (Drag force)
***************************************************************************/
#include "../include/cannon_aero.h"
#include "trick_utils/math/include/trick_math.h"
int cannon_force_drag(
CANNON_AERO *C )
{
double k ;
double speed ;
speed = V_MAG( C->vel ) ;
/* k = -1/2*rho*Cd*A*|V| */
k = (-0.5)*C->air_density*C->coefficient_drag*C->ball_area*speed ;
/* Force_drag = k*V = -1/2*rho*Cd*A*|V|*V */
V_SCALE( C->force_drag, C->vel, k ) ;
C->mag_force_drag = V_MAG( C->force_drag ) ;
return 0 ;
}

15
trick_models/baseball/aero/src/cannon_force_gravity.c
View File

@ -1,15 +0,0 @@
/*********************************** TRICK HEADER **************************
PURPOSE: ( Gravitational force on cannonball )
***************************************************************************/
#include "../include/cannon_aero.h"
int cannon_force_gravity(
CANNON_AERO *C )
{
C->force_gravity[0] = 0.0 ;
C->force_gravity[1] = 0.0 ;
C->force_gravity[2] = C->mass*C->g ;
return 0 ;
}

18
trick_models/baseball/aero/src/cannon_force_jet.c
View File

@ -1,18 +0,0 @@
/*********************************** TRICK HEADER **************************
PURPOSE: ( Jet fire force )
***************************************************************************/
#include "../include/cannon_aero.h"
int cannon_force_jet(
CANNON_AERO *C )
{
if ( C->jet_on && C->jet_count < 4 ) {
C->force_jet[2] = C->force_jet_Z_plus ;
C->jet_count++ ;
C->jet_on = 0 ;
} else {
C->force_jet[2] = 0.0 ;
}
return 0 ;
}

62
trick_models/baseball/aero/src/cannon_force_lift.c
View File

@ -1,62 +0,0 @@
/*********************************** TRICK HEADER **************************
PURPOSE: (Lift-Magnus Force)
***************************************************************************/
#include "../include/cannon_aero.h"
#include "trick_utils/math/include/trick_math.h"
int cannon_force_lift(
CANNON_AERO *C )
{
double w_cross_v[3] ; double norm_w_cross_v[3] ;
double k, speed ;
speed = V_MAG( C->vel ) ;
if ( speed != 0.0 ) {
C->spin_parameter = C->ball_radius*V_MAG( C->omega )/speed ;
} else {
C->spin_parameter = 0.0000000001 ;
}
if ( C->lift_method == Smits_Smith ) {
C->coefficient_lift = (0.54)*pow(C->spin_parameter, 0.4) ;
}
if ( C->lift_method == Tombras ) {
C->coefficient_lift = 1/( 2.022 + 0.981*speed/V_MAG( C->omega));
}
switch ( C->lift_method ) {
case Hard_Coded_Coefficient_Lift:
case Smits_Smith:
case Tombras:
/* k = 1/2*rho*Cl*A*V^2 */
k = (0.5)*C->air_density*C->coefficient_lift*
C->ball_area*speed*speed ;
/* F = k*(w x V)/|w x V| */
V_CROSS( w_cross_v, C->omega, C->vel ) ;
V_NORM( norm_w_cross_v, w_cross_v ) ;
V_SCALE( C->force_magnus, norm_w_cross_v, k) ;
break ;
case Adair_Giordano:
V_CROSS( w_cross_v, C->omega, C->vel ) ;
k = 0.00041*C->mass ;
V_SCALE( C->force_magnus, w_cross_v, k) ;
/* Backwards calculation for Cl */
C->coefficient_lift = (2*k*V_MAG(w_cross_v))/
(C->air_density*C->ball_area*speed*speed) ;
break ;
}
C->mag_force_magnus = V_MAG( C->force_magnus ) ;
return 0 ;
}

47
trick_models/baseball/aero/src/cannon_impact_aero.c
View File

@ -1,47 +0,0 @@
/*********************************** TRICK HEADER **************************
PURPOSE: (Kaboom!!!)
***************************************************************************/
#include <stdio.h>
#include "../include/cannon_aero.h"
#include "sim_services/include/executive.h"
#include "sim_services/include/exec_proto.h"
#include "sim_services/include/regula_falsi.h"
#include "sim_services/include/dr_proto.h"
double cannon_impact_aero(
CANNON_AERO* C,
double* time,
int *event_evaluate_tgo )
{
double tgo ;
EXECUTIVE* E ;
E = exec_get_exec();
if( *event_evaluate_tgo ) {
/* Calculate time to go before impact */
C->rf.error = C->pos[0] ;
tgo = regula_falsi( *time , &(C->rf) ) ;
} else {
/* Ball impact */
reset_regula_falsi( *time , &(C->rf) ) ;
tgo = 0.0 ;
C->vel[0] = 0.0 ; C->vel[1] = 0.0 ; C->vel[2] = 0.0 ;
C->acc[0] = 0.0 ; C->acc[1] = 0.0 ; C->acc[2] = 0.0 ;
C->g = 0.0 ;
fprintf(stderr, "Impact time, pos : %.9lf %.9lf\n",
*time, C->pos[0] );
dr_record_binary( &E->record.group[0], *time);
}
return( tgo ) ;
}

47
trick_models/baseball/aero/src/cannon_impact_monte.c
View File

@ -1,47 +0,0 @@
/*********************************** TRICK HEADER **************************
PURPOSE: (Kaboom!!!)
***************************************************************************/
#include <stdio.h>
#include "../include/cannon_aero.h"
#include "sim_services/include/executive.h"
#include "sim_services/include/exec_proto.h"
#include "sim_services/include/regula_falsi.h"
#include "sim_services/include/dr_proto.h"
double cannon_impact_monte(
CANNON_AERO* C,
double* time,
int *event_evaluate_tgo )
{
double tgo ;
EXECUTIVE* E ;
E = exec_get_exec();
if( *event_evaluate_tgo ) {
/* Calculate time to go before impact */
C->rf.error = C->pos[2] ;
tgo = regula_falsi( *time , &(C->rf) ) ;
} else {
/* Ball impact */
reset_regula_falsi( *time , &(C->rf) ) ;
tgo = 0.0 ;
C->vel[0] = 0.0 ; C->vel[1] = 0.0 ; C->vel[2] = 0.0 ;
C->acc[0] = 0.0 ; C->acc[1] = 0.0 ; C->acc[2] = 0.0 ;
C->g = 0.0 ;
fprintf(stderr, "Impact time, pos : %.9lf %.9lf\n",
*time, C->pos[0] );
dr_record_binary( &E->record.group[0], *time);
}
return( tgo ) ;
}

19
trick_models/baseball/aero/src/cannon_init_aero.c
View File

@ -1,19 +0,0 @@
/*********************************** TRICK HEADER **************************
PURPOSE: ( Cannon initialization )
***************************************************************************/
#include <math.h>
#include <stdio.h>
#include "../include/cannon_aero.h"
#include "trick_utils/math/include/trick_math.h"
int cannon_init_aero(
CANNON_AERO* C)
{
/* Convert omega from spherical (almost) to rectangular */
C->omega[0] = C->omega0*sin(M_PI/2.0 - C->phi)*cos(C->theta) ;
C->omega[1] = C->omega0*sin(M_PI/2.0 - C->phi)*sin(C->theta) ;
C->omega[2] = C->omega0*cos(M_PI/2.0 - C->phi) ;
return ( 0 );
}

55
trick_models/baseball/aero/src/cannon_integ_aero.c
View File

@ -1,55 +0,0 @@
/*********************************** TRICK HEADER **************************
PURPOSE: (Cannon integration)
***************************************************************************/
#include "sim_services/Integrator/include/integrator_c_intf.h"
#include "../include/cannon_aero.h"
int cannon_integ_aero(
CANNON_AERO* C)
{
int ipass;
/* GET SHORTHAND NOTATION FOR DATA STRUCTURES */
CANNON_AERO_OUT *CAO = &(C->out) ;
/* LOAD THE POSITION AND VELOCITY STATES */
load_state(
&CAO->position[0] ,
&CAO->position[1] ,
&CAO->position[2] ,
&CAO->velocity[0] ,
&CAO->velocity[1] ,
&CAO->velocity[2] ,
NULL
);
/* LOAD THE POSITION AND VELOCITY STATE DERIVATIVES */
load_deriv(
&CAO->velocity[0] ,
&CAO->velocity[1] ,
&CAO->velocity[2] ,
&CAO->acceleration[0] ,
&CAO->acceleration[1] ,
&CAO->acceleration[2] ,
NULL
);
/* CALL THE TRICK INTEGRATION SERVICE */
ipass = integrate();
/* UNLOAD THE NEW POSITION AND VELOCITY STATES */
unload_state(
&CAO->position[0] ,
&CAO->position[1] ,
&CAO->position[2] ,
&CAO->velocity[0] ,
&CAO->velocity[1] ,
&CAO->velocity[2] ,
NULL
);
/* RETURN */
return( ipass );
}

25
trick_models/baseball/aero/src/cannon_jet_control.c
View File

@ -1,25 +0,0 @@
/*********************************** TRICK HEADER **************************
PURPOSE: (Controls when jets are fired)
***************************************************************************/
#include "../include/cannon_aero.h"
#include "sim_services/include/exec_proto.h"
#include "trick_utils/math/include/trick_math.h"
int cannon_jet_control( CANNON_AERO* C )
{
double sim_time ;
#define CANNON_EQUALS(X,Y) ( fabs(X - Y) < 1.0e-9 ) ? 1 : 0
sim_time = exec_get_sim_time();
if ( CANNON_EQUALS(sim_time, roundoff(0.1, C->time_to_fire_jet_1)) ||
CANNON_EQUALS(sim_time, roundoff(0.1, C->time_to_fire_jet_2)) ||
CANNON_EQUALS(sim_time, roundoff(0.1, C->time_to_fire_jet_3)) ||
CANNON_EQUALS(sim_time, roundoff(0.1, C->time_to_fire_jet_4)) ) {
C->jet_on = 1 ;
}
return(0) ;
}

67
trick_models/baseball/optim/include/amoeba.h
View File

@ -1,67 +0,0 @@
/*
* Vetter's implementation of the Nelder-Mead "amoeba" algorithm.
* Was made strictly for Trick tutorial. Was a whim.
* Reference for algorithm (not code):
* http://www.research.ibm.com/infoecon/paps/html/amec99_bundle/node8.html
*/
#ifndef AMOEBA_H
#define AMOEBA_H
#define AMOEBA_ALPHA 1.0 /* reflection constant */
#define AMOEBA_BETA 1.0 /* expansion constant */
#define AMOEBA_ZETA 0.5 /* contraction constant */
#define AMOEBA_ETA 0.5 /* shrinkage constant */
typedef enum {
VERTICES,
CALC_CENTROID_POINT,
CALC_REFLECTION_POINT,
REFLECT,
EXPAND,
CONTRACT,
SHRINK,
} AMOEBA_STATE ;
typedef struct {
int debug ; /* -- Turn on some printing */
int max_steps ; /* -- Even if we don't find solution, break */
int num_dims ; /* -- Num dims of each vertice */
int num_vertices ; /* -- Num vertices in simplex */
double epsilon ; /* -- Bail if max dist simplex < epsilon */
double** vertices ; /* -- Rows are vectors of num_dims */
double** init_simplex ; /* -- Initial simplex */
double* x_cent ; /* -- Center of vertices-exclude worst vertex */
double* x_refl ; /* -- Point of reflection */
double* x_expa ; /* -- Point of expansion */
double* x_cont ; /* -- Point of contraction */
double* curr_point ; /* -- Current point of interest for sim use:
Simplex vertice, centroid, reflection,
expansion or contraction point */
int curr_vertex ; /* -- Current vertex for sim state machine */
AMOEBA_STATE state ; /* -- For sim's amoeba state machine */
} AMOEBA ;
/* Prototypes */
void amoeba_init( AMOEBA* A, int num_dims, double epsilon, int max_steps,
double* simplex_point, double simplex_size ) ;
void amoeba_print( AMOEBA* A ) ;
void amoeba_print_point( int num_dims, double* point ) ;
void amoeba_go( AMOEBA* A ) ;
double my_func( double* point ) ;
void amoeba_quit( AMOEBA* A ) ;
void amoeba_order( AMOEBA* A ) ;
int amoeba_reflect( AMOEBA* A ) ;
void amoeba_calc_reflection_point( AMOEBA* A ) ;
int amoeba_expand( AMOEBA* A ) ;
void amoeba_calc_expansion_point( AMOEBA* A ) ;
int amoeba_contract( AMOEBA* A ) ;
void amoeba_calc_contraction_point( AMOEBA* A ) ;
void amoeba_shrink( AMOEBA* A ) ;
int amoeba_satisfied( AMOEBA* A ) ;
double amoeba_distance_between( int num_dims, double* x, double* y) ;
int amoeba_compare_vertices( const void* a, const void* b) ;
#endif

566
trick_models/baseball/optim/src/amoeba.c
View File

@ -1,566 +0,0 @@
/*********************************** TRICK HEADER **************************
PURPOSE: ( Amoeba )
***************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "../include/amoeba.h"
double my_func( double* point ) {
/* In the case of the simulation, the function value is
* actually stored in the point
*/
int num_dims = 4 ;
return( point[num_dims] ) ;
}
void amoeba_go( AMOEBA* A ) {
int step ;
step = 1 ;
while( ! amoeba_satisfied(A) && step < A->max_steps ) {
step++ ;
amoeba_order( A ) ;
if ( amoeba_reflect( A ) ) {
continue ;
}
if ( amoeba_expand( A ) ) {
continue ;
}
if ( amoeba_contract( A ) ) {
continue ;
}
amoeba_shrink( A ) ;
}
amoeba_order( A ) ;
fprintf(stderr,"\n\nAmoeba search concluded.\n") ;
fprintf(stderr,"Number iterations: %d \n", step) ;
fprintf(stderr,"x=(");
amoeba_print_point( A->num_dims, A->vertices[0] ) ;
fprintf(stderr,")\n");
fprintf(stderr,"F(x)=%.6lf\n\n", my_func(A->vertices[0]));
}
void amoeba_shrink( AMOEBA* A ) {
int ii, jj ;
double* x ;
double* x_0 ;
double* x_n ;
x_0 = A->vertices[0] ;
x_n = A->vertices[A->num_dims] ;
if ( A->debug ) {
fprintf(stderr,"Shrinking simplex...\n");
amoeba_print( A ) ;
fprintf(stderr,"\n\n");
}
for ( ii = 0 ; ii < A->num_vertices ; ii++ ) {
x = A->vertices[ii] ;
for ( jj = 0 ; jj < A->num_dims ; jj++ ) {
x[jj] = x_0[jj] + AMOEBA_ETA*( x[jj] - x_0[jj] ) ;
}
}
if ( A->debug ) {
fprintf(stderr,"Shrank simplex.\n");
amoeba_print( A ) ;
fprintf(stderr,"\n\n");
}
}
void amoeba_calc_contraction_point( AMOEBA* A ) {
int ii ;
double* x_c ;
double* x_n ;
double* x_cont ;
/* Short hand */
x_c = A->x_cent ;
x_n = A->vertices[A->num_vertices-1] ;
x_cont = A->x_cont ;
/* Calculate contraction point
* x_contract = x_center + zeta*(x_center - x_worst)
*/
for ( ii = 0 ; ii < A->num_dims ; ii++ ) {
x_cont[ii] = x_c[ii] + AMOEBA_ZETA*(x_c[ii] - x_n[ii]) ;
}
}
int amoeba_contract( AMOEBA* A ) {
int ii ;
double* x_r ;
double* x_n ;
double* x_n_1 ;
double* x_c ;
double* x_cont ;
double* x_e ;
/* Short hand */
x_c = A->x_cent ;
x_r = A->x_refl ;
x_n = A->vertices[A->num_vertices-1] ;
x_n_1 = A->vertices[A->num_vertices-2] ;
x_cont = A->x_cont ;
x_e = A->x_expa ;
if ( A->debug ) {
fprintf(stderr,">> Try contraction.\n");
}
/* Contract if F(x_r) < F(x_n-1) */
if ( my_func(x_r) >= my_func(x_n_1) ) {
if ( A->debug ) {
fprintf(stderr,">> Reject contraction.\n");
}
return 0 ;
}
/* Calculate contraction point
* x_contract = x_center + zeta*(x_center - x_worst)
*/
for ( ii = 0 ; ii < A->num_dims ; ii++ ) {
x_cont[ii] = x_c[ii] + AMOEBA_ZETA*(x_c[ii] - x_n[ii]) ;
}
if ( my_func(x_cont) >= my_func(x_n) ) {
/* Accept contraction point */
if ( A->debug ) {
fprintf(stderr,">> Accept contraction about x_cont(");
amoeba_print_point( A->num_dims, x_cont ) ;
fprintf(stderr,") F(x_cont)=%.2lf \n",
my_func(x_cont) ) ;
}
for ( ii = 0 ; ii < A->num_dims+1 ; ii++ ) {
x_n[ii] = x_cont[ii] ;
}
return 1 ;
} else {
/* Reject contraction point */
if ( A->debug ) {
fprintf(stderr,">> Reject contraction about x_cont(");
amoeba_print_point( A->num_dims, x_cont ) ;
fprintf(stderr,") F(x_cont)=%.2lf \n",
my_func(x_cont) ) ;
}
return 0 ;
}
}
void amoeba_calc_expansion_point( AMOEBA* A ) {
int ii ;
double* x_r ;
double* x_c ;
double* x_e ;
/* Short hand */
x_r = A->x_refl ;
x_c = A->x_cent ;
x_e = A->x_expa ;
/* Calculate expansion point
* x_expa = x_refl + beta*(x_refl - x_cent)
*/
for ( ii = 0 ; ii < A->num_dims ; ii++ ) {
x_e[ii] = x_r[ii] + AMOEBA_BETA*(x_r[ii] - x_c[ii]) ;
}
}
int amoeba_expand( AMOEBA* A ) {
int ii ;
double* x_0 ;
double* x_r ;
double* x_n ;
double* x_c ;
double* x_e ;
/* Short hand */
x_0 = A->vertices[0] ;
x_r = A->x_refl ;
x_n = A->vertices[A->num_vertices-1] ;
x_c = A->x_cent ;
x_e = A->x_expa ;
if ( A->debug ) {
fprintf(stderr,">> Try expansion.\n");
}
/* Expand if F(x_r) > F(x_0) */
if ( my_func(x_r) <= my_func(x_0) ) {
if ( A->debug ) {
fprintf(stderr,">> Reject expansion.\n");
}
return 0 ;
}
if ( A->debug ) {
fprintf(stderr,">> Accept expansion.\n");
}
if ( my_func(x_e) >= my_func(x_r) ) {
/* Accept expansion point */
if ( A->debug ) {
fprintf(stderr,"Accept x_e(");
amoeba_print_point( A->num_dims, x_e ) ;
fprintf(stderr,") F(x_e)=%.2lf \n", my_func(x_e) ) ;
}
for ( ii = 0 ; ii < A->num_dims+1 ; ii++ ) {
x_n[ii] = x_e[ii] ;
}
} else {
/* Accept reflection point */
if ( A->debug ) {
fprintf(stderr,"Accept x_r(");
amoeba_print_point( A->num_dims, x_r ) ;
fprintf(stderr,") F(x_r)=%.2lf \n", my_func(x_r) ) ;
}
for ( ii = 0 ; ii < A->num_dims+1 ; ii++ ) {
x_n[ii] = x_r[ii] ;
}
}
return 1 ;
}
double amoeba_distance_between( int num_dims, double* x, double* y ) {
int ii ;
double sum_squares ;
sum_squares = 0 ;
for ( ii = 0 ; ii < num_dims ; ii++ ) {
sum_squares += (x[ii] - y[ii])*(x[ii] - y[ii]);
}
return ( sqrt(sum_squares) ) ;
}
int amoeba_satisfied( AMOEBA* A ) {
/*
* Find max distance between points in the simplex
* If under epsilon, amoeba satisfied
* If this isn't satisfied, see if F(x) satisfied within epsilon
*/
int ii, jj ;
double dist, max_dist ;
double* x ;
double* y ;
double f1 ;
double f2 ;
/*
* Check max dist between points in simplex
* If the amoeba gets too small, no point in shrinking
* infintessimally (sp?).
*/
max_dist = 0.0 ;
for ( ii = 0 ; ii < A->num_vertices - 1 ; ii++ ) {
x = A->vertices[ii] ;
for ( jj = ii+1 ; jj < A->num_vertices ; jj++ ) {
y = A->vertices[jj] ;
dist = amoeba_distance_between( A->num_dims, x, y ) ;
if ( dist > max_dist ) {
max_dist = dist ;
}
}
}
if ( max_dist < A->epsilon ) {
return 1 ;
}
/*
* See if diff between function vals less than epsilon
*/
max_dist = 0 ;
for ( ii = 0 ; ii < A->num_vertices - 1 ; ii++ ) {
f1 = A->vertices[ii][A->num_dims] ;
for ( jj = ii+1 ; jj < A->num_vertices ; jj++ ) {
f2 = A->vertices[jj][A->num_dims] ;
dist = fabs( f1 - f2 ) ;
if ( dist > max_dist ) {
max_dist = dist ;
}
}
}
if ( max_dist < A->epsilon ) {
return 1 ;
}
/*
* Unsatisfied amoeba :-(
*/
return 0 ;
}
void amoeba_calc_reflection_point( AMOEBA* A ) {
int ii ;
double* x_r ;
double* x_n ;
double* x_c ;
/* Short hand */
x_r = A->x_refl ;
x_n = A->vertices[A->num_vertices-1] ;
x_c = A->x_cent ;
/* Calculate reflection point
* x_refl = x_cent + alpha*(x_cent - x_worst)
*/
for ( ii = 0 ; ii < A->num_dims ; ii++ ) {
x_r[ii] = x_c[ii] + AMOEBA_ALPHA*(x_c[ii] - x_n[ii]) ;
}
}
int amoeba_reflect( AMOEBA* A ) {
int ii ;
double* x_0 ;
double* x_r ;
double* x_n ;
double* x_c ;
if ( A->debug ) {
fprintf(stderr,">> Try reflection.\n");
}
/* Short hand */
x_0 = A->vertices[0] ;
x_r = A->x_refl ;
x_n = A->vertices[A->num_vertices-1] ;
x_c = A->x_cent ;
if ( my_func(x_0) >= my_func(x_r) && my_func(x_r) > my_func(x_n) ) {
/* Accept reflection point --- replace worst point with x_r */
if ( A->debug ) {
fprintf(stderr,"Accept reflection x_r(");
amoeba_print_point( A->num_dims, x_r ) ;
fprintf(stderr,") F(x_r)=%.2lf \n", my_func(x_r) ) ;
}
for ( ii = 0 ; ii < A->num_dims+1 ; ii++ ) {
x_n[ii] = x_r[ii] ;
}
return 1 ;
} else {
if ( A->debug ) {
fprintf(stderr,"Reject reflection x_r(");
amoeba_print_point( A->num_dims, x_r ) ;
fprintf(stderr,") F(x_r)=%.2lf \n", my_func(x_r) ) ;
}
return 0 ;
}
}
void amoeba_order( AMOEBA* A ) {
int ii, jj ;
/* Order vertices based on function results.
* Order from best to worst (i.e. max to min)
*/
qsort( A->vertices,
A->num_vertices, sizeof(double*),
amoeba_compare_vertices) ;
/* Calculate centroid --- excluding worst point */
for ( jj = 0 ; jj < A->num_dims ; jj++ ) {
A->x_cent[jj] = 0 ;
}
for ( jj = 0 ; jj < A->num_dims ; jj++ ) {
for ( ii = 0 ; ii < A->num_vertices - 1 ; ii++ ) {
A->x_cent[jj] += A->vertices[ii][jj] ;
}
A->x_cent[jj] /= (double)(A->num_vertices - 1.0 ) ;
}
if ( A->debug ) {
fprintf(stderr,"Ordered Amoeba.\n");
amoeba_print( A ) ;
fprintf(stderr,"\n");
}
}
int amoeba_compare_vertices( const void* ap, const void* bp ) {
double a, b ;
a = my_func( *((double**) ap) ) ;
b = my_func( *((double**) bp) ) ;
if ( a > b ) {
return -1 ;
} else if ( a < b ) {
return 1 ;
} else {
return 0 ;
}
}
/*
* num_dims: Number of dimensions in domain space.
* Example: If the function we are maximizing is F(x,y,z)
* The number of dims is 3.
*
* epsilon: The error tolerance before amoeba algorithm stops
* A good number is something like: 1.0e-9
*
* max_steps: Even if we don't reach a solution by "max_steps", stop.
*
* We have to create a first simplex to begin amoeba algorithm.
* There are two ways to do this.
* First : Send simplex_point and simplex_size
* This defines a simplex with one vertex at "simplex_point"
* All other points are spaced "simplex_size" in the unit vec direction
* Second : If "simplex_point" is NULL, it is assumed user setup
* A->init_simplex before entering amoeba_init.
*/
void amoeba_init( AMOEBA* A, int num_dims, double epsilon, int max_steps,
double* simplex_point, double simplex_size ) {
int ii, jj ;
if ( A->init_simplex == NULL && simplex_point == NULL ) {
fprintf(stderr, "Initial simplex undefined by user.\n"
"Set A->init_simplex or "
"simplex_point and simplex_size\n") ;
exit(-1) ;
}
A->state = VERTICES ;
A->num_dims = num_dims ;
A->num_vertices = num_dims + 1 ;
A->curr_vertex = 0 ;
A->epsilon = epsilon ;
A->max_steps = max_steps ;
A->vertices = (double**) malloc (A->num_vertices*sizeof(double));
for ( ii = 0 ; ii < A->num_vertices ; ii++ ) {
A->vertices[ii] = (double*) malloc
((A->num_dims+1)*sizeof(double));
}
if ( simplex_point != NULL ) {
/*
* Use "simplex_point" as one vertice of first simplex. Then
* move in "simplex_size" along unit vectors to create the
* rest of the first guess simplex.
*
* First_simplex = simplex_point + simplex_size*unit_vecs
*/
for ( jj = 0 ; jj < A->num_dims ; jj++ ) {
A->vertices[0][jj] = simplex_point[jj] ;
}
/* Load function vals with zero */
A->vertices[0][jj] = 0.0 ;
for ( ii = 1 ; ii < A->num_vertices ; ii++ ) {
for ( jj = 0 ; jj < A->num_dims ; jj++ ) {
if ( jj+1 == ii ) {
A->vertices[ii][jj] =
simplex_point[jj] - simplex_size ;
} else {
A->vertices[ii][jj] = simplex_point[jj] ;
}
}
/* Load function vals with zero */
A->vertices[ii][jj] = 0.0 ;
}
} else {
/* Assume user has defined first simplex */
for ( ii = 0 ; ii < A->num_vertices ; ii++ ) {
for ( jj = 0 ; jj < A->num_dims ; jj++ ) {
A->vertices[ii][jj] = A->init_simplex[ii][jj] ;
}
/* Load function vals with zero */
A->vertices[ii][jj] = 0.0 ;
}
}
A->x_cent = (double*) malloc ( (A->num_dims+1) * sizeof(double) ) ;
A->x_refl = (double*) malloc ( (A->num_dims+1) * sizeof(double) ) ;
A->x_expa = (double*) malloc ( (A->num_dims+1) * sizeof(double) ) ;
A->x_cont = (double*) malloc ( (A->num_dims+1) * sizeof(double) ) ;
/* Init func results */
A->x_cent[A->num_dims] = 0.0 ;
A->x_refl[A->num_dims] = 0.0 ;
A->x_expa[A->num_dims] = 0.0 ;
A->x_cont[A->num_dims] = 0.0 ;
if ( A->debug ) {
fprintf(stderr,"Initial simplex with function evals.\n");
amoeba_print( A ) ;
fprintf(stderr,"\n\n");
}
}
void amoeba_quit( AMOEBA* A ) {
int ii ;
for ( ii = 0 ; ii < A->num_vertices ; ii++ ) {
free( A->vertices[ii] ) ;
}
free( A->vertices ) ;
free( A->x_cent ) ;
free( A->x_refl ) ;
free( A->x_expa ) ;
free( A->x_cont ) ;
}
void amoeba_print_point( int num_dims, double* point ) {
int ii;
for ( ii = 0 ; ii < num_dims-1 ; ii++ ) {
fprintf(stderr,"%.2lf ", point[ii]);
}
fprintf(stderr,"%.2lf", point[num_dims-1]);
}
void amoeba_print( AMOEBA* A ) {
int ii, jj ;
for ( ii = 0 ; ii < A->num_vertices ; ii++ ) {
for ( jj = 0 ; jj < A->num_dims ; jj++ ) {
fprintf(stderr,"%.2lf ", A->vertices[ii][jj] ) ;
}
fprintf(stderr," > %.2lf \n", A->vertices[ii][jj]);
}
}

15
trick_models/baseball/optim/src/cannon_init_amoeba.c
View File

@ -1,15 +0,0 @@
/*********************************** TRICK HEADER **************************
PURPOSE: ( Init the Amoeba )
LIBRARY DEPENDENCY: ((amoeba.o))
PROGRAMMERS: (Keith)
***************************************************************************/
#include <stdio.h>
#include "cannon/aero/include/cannon_aero.h"
#include "../include/amoeba.h"
int cannon_init_amoeba(
AMOEBA* A )
{
amoeba_init( A, 4, 1.0e-3, 100, NULL, 0.0 ) ;
return ( 0 );
}

35
trick_models/baseball/optim/src/cannon_post_master.c
View File

@ -1,35 +0,0 @@
/*********************************** TRICK HEADER **************************
PURPOSE: (Read slave's sim evaluation)
LIBRARY DEPENDENCY: ((amoeba.o))
PROGRAMMER: ((keith))
***************************************************************************/
#include <stdio.h>
#include "cannon/aero/include/cannon_aero.h"
#include "../include/amoeba.h"
#include "sim_services/include/executive.h"
#include "sim_services/include/exec_proto.h"
int cannon_post_master(
CANNON_AERO* C,
AMOEBA* A )
{
CANNON_AERO C_curr ;
EXECUTIVE* E ;
E = exec_get_exec();
/* Read slave's results */
tc_read( &E->monte.work.data_conn,
(char*) &C_curr, sizeof(CANNON_AERO) );
fprintf(stderr, "%03d> F(", E->monte.work.curr_run);
amoeba_print_point(4, A->curr_point) ;
fprintf(stderr, ") = %.6lf\n", C_curr.pos[0]);
/*
* Load function result for either:
* simplex vertice, centroid, reflection, contraction or expansion point
*/
A->curr_point[A->num_dims] = C_curr.pos[0] ;
return(0) ;
}

20
trick_models/baseball/optim/src/cannon_post_slave.c
View File

@ -1,20 +0,0 @@
/*********************************** TRICK HEADER **************************
PURPOSE: (Get slave sim's evaluation of x)
***************************************************************************/
#include "cannon/aero/include/cannon_aero.h"
#include "sim_services/include/executive.h"
#include "sim_services/include/exec_proto.h"
int cannon_post_slave(
/* RETURN: -- Always return zero */
CANNON_AERO* C) /* INOUT: -- Parameter */
{
EXECUTIVE* E ;
E = exec_get_exec();
/* Send slave results */
tc_write( &E->monte.work.data_conn, (char*) C, sizeof(CANNON_AERO) );
return(0) ;
}

105
trick_models/baseball/optim/src/cannon_pre_master.c
View File

@ -1,105 +0,0 @@
/*********************************** TRICK HEADER **************************
PURPOSE: (Master optimization)
LIBRARY_DEPENDENCY: (amoeba.o)
PROGRAMMER: (keith)
***************************************************************************/
#include <stdio.h>
#include "cannon/aero/include/cannon_aero.h"
#include "../include/amoeba.h"
#include "sim_services/include/executive.h"
#include "sim_services/include/exec_proto.h"
int cannon_pre_master(
/* RETURN: -- Always return zero */
CANNON_AERO* C, /* INOUT: -- Parameter */
AMOEBA* A)
{
while ( 1 ) {
switch ( A->state ) {
case VERTICES:
A->curr_point = A->vertices[A->curr_vertex] ;
if ( A->curr_vertex == A->num_vertices ) {
A->state = CALC_CENTROID_POINT ;
} else {
fprintf(stderr, "V[%d] ", A->curr_vertex);
}
A->curr_vertex++ ;
break ;
case CALC_CENTROID_POINT:
fprintf(stderr, "CENT ");
amoeba_order( A ) ;
A->curr_point = A->x_cent ;
A->state = CALC_REFLECTION_POINT ;
break ;
case CALC_REFLECTION_POINT:
fprintf(stderr, "REFL ");
amoeba_calc_reflection_point( A ) ;
A->curr_point = A->x_refl ;
A->state = REFLECT ;
break ;
case REFLECT:
if ( amoeba_reflect( A ) ) {
A->state = CALC_CENTROID_POINT ;
} else {
fprintf(stderr, "EXPA ");
amoeba_calc_expansion_point( A ) ;
A->curr_point = A->x_expa ;
A->state = EXPAND ;
}
break ;
case EXPAND:
if ( amoeba_expand( A ) ) {
A->state = CALC_CENTROID_POINT ;
} else {
fprintf(stderr, "CONT ");
amoeba_calc_contraction_point( A ) ;
A->curr_point = A->x_cont ;
A->state = CONTRACT ;
}
break ;
case CONTRACT:
if ( amoeba_contract( A ) ) {
A->state = CALC_CENTROID_POINT ;
} else {
A->state = SHRINK ;
}
break ;
case SHRINK:
amoeba_shrink( A ) ;
fprintf(stderr, "V[0] ");
A->curr_point = A->vertices[0] ;
A->curr_vertex = 1 ;
A->state = VERTICES ;
break ;
}
if ( amoeba_satisfied( A ) && A->state != VERTICES ) {
exec_terminate( "cannon_pre_master",
"Amoeba has found a solution." ) ;
}
if ( A->state == CALC_CENTROID_POINT ||
A->state == SHRINK ) {
continue ;
} else {
break ;
}
}
C->time_to_fire_jet_1 = A->curr_point[0] ;
C->time_to_fire_jet_2 = A->curr_point[1] ;
C->time_to_fire_jet_3 = A->curr_point[2] ;
C->time_to_fire_jet_4 = A->curr_point[3] ;
return(0) ;
}

23
trick_models/baseball/optim/src/cannon_pre_slave.c
View File

@ -1,23 +0,0 @@
/*********************************** TRICK HEADER **************************
PURPOSE: (Slave optimization)
LIBRARY_DEPENDENCY: (amoeba.o)
PROGRAMMER: (keith)
***************************************************************************/
#include "cannon/aero/include/cannon_aero.h"
#include "sim_services/include/executive.h"
#include "sim_services/include/exec_proto.h"
#include <stdio.h>
int cannon_pre_slave(
/* RETURN: -- Always return zero */
CANNON_AERO* C) /* INOUT: -- Parameter */
{
fprintf(stderr, "SLAVE POINT: (%.6lf %.6lf %.6lf %.6lf)\n",
C->time_to_fire_jet_1,
C->time_to_fire_jet_2,
C->time_to_fire_jet_3,
C->time_to_fire_jet_4) ;
return(0) ;
}