001/* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017 018package org.apache.commons.math3.ode.events; 019 020import java.util.Arrays; 021 022/** Wrapper used to detect only increasing or decreasing events. 023 * 024 * <p>General {@link EventHandler events} are defined implicitly 025 * by a {@link EventHandler#g(double, double[]) g function} crossing 026 * zero. This function needs to be continuous in the event neighborhood, 027 * and its sign must remain consistent between events. This implies that 028 * during an ODE integration, events triggered are alternately events 029 * for which the function increases from negative to positive values, 030 * and events for which the function decreases from positive to 031 * negative values. 032 * </p> 033 * 034 * <p>Sometimes, users are only interested in one type of event (say 035 * increasing events for example) and not in the other type. In these 036 * cases, looking precisely for all events location and triggering 037 * events that will later be ignored is a waste of computing time.</p> 038 * 039 * <p>Users can wrap a regular {@link EventHandler event handler} in 040 * an instance of this class and provide this wrapping instance to 041 * the {@link org.apache.commons.math3.ode.FirstOrderIntegrator ODE solver} 042 * in order to avoid wasting time looking for uninteresting events. 043 * The wrapper will intercept the calls to the {@link 044 * EventHandler#g(double, double[]) g function} and to the {@link 045 * EventHandler#eventOccurred(double, double[], boolean) 046 * eventOccurred} method in order to ignore uninteresting events. The 047 * wrapped regular {@link EventHandler event handler} will the see only 048 * the interesting events, i.e. either only {@code increasing} events or 049 * {@code decreasing} events. the number of calls to the {@link 050 * EventHandler#g(double, double[]) g function} will also be reduced.</p> 051 * 052 * @since 3.2 053 */ 054 055public class EventFilter implements EventHandler { 056 057 /** Number of past transformers updates stored. */ 058 private static final int HISTORY_SIZE = 100; 059 060 /** Wrapped event handler. */ 061 private final EventHandler rawHandler; 062 063 /** Filter to use. */ 064 private final FilterType filter; 065 066 /** Transformers of the g function. */ 067 private final Transformer[] transformers; 068 069 /** Update time of the transformers. */ 070 private final double[] updates; 071 072 /** Indicator for forward integration. */ 073 private boolean forward; 074 075 /** Extreme time encountered so far. */ 076 private double extremeT; 077 078 /** Wrap an {@link EventHandler event handler}. 079 * @param rawHandler event handler to wrap 080 * @param filter filter to use 081 */ 082 public EventFilter(final EventHandler rawHandler, final FilterType filter) { 083 this.rawHandler = rawHandler; 084 this.filter = filter; 085 this.transformers = new Transformer[HISTORY_SIZE]; 086 this.updates = new double[HISTORY_SIZE]; 087 } 088 089 /** {@inheritDoc} */ 090 public void init(double t0, double[] y0, double t) { 091 092 // delegate to raw handler 093 rawHandler.init(t0, y0, t); 094 095 // initialize events triggering logic 096 forward = t >= t0; 097 extremeT = forward ? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY; 098 Arrays.fill(transformers, Transformer.UNINITIALIZED); 099 Arrays.fill(updates, extremeT); 100 101 } 102 103 /** {@inheritDoc} */ 104 public double g(double t, double[] y) { 105 106 final double rawG = rawHandler.g(t, y); 107 108 // search which transformer should be applied to g 109 if (forward) { 110 final int last = transformers.length - 1; 111 if (extremeT < t) { 112 // we are at the forward end of the history 113 114 // check if a new rough root has been crossed 115 final Transformer previous = transformers[last]; 116 final Transformer next = filter.selectTransformer(previous, rawG, forward); 117 if (next != previous) { 118 // there is a root somewhere between extremeT and t. 119 // the new transformer is valid for t (this is how we have just computed 120 // it above), but it is in fact valid on both sides of the root, so 121 // it was already valid before t and even up to previous time. We store 122 // the switch at extremeT for safety, to ensure the previous transformer 123 // is not applied too close of the root 124 System.arraycopy(updates, 1, updates, 0, last); 125 System.arraycopy(transformers, 1, transformers, 0, last); 126 updates[last] = extremeT; 127 transformers[last] = next; 128 } 129 130 extremeT = t; 131 132 // apply the transform 133 return next.transformed(rawG); 134 135 } else { 136 // we are in the middle of the history 137 138 // select the transformer 139 for (int i = last; i > 0; --i) { 140 if (updates[i] <= t) { 141 // apply the transform 142 return transformers[i].transformed(rawG); 143 } 144 } 145 146 return transformers[0].transformed(rawG); 147 148 } 149 } else { 150 if (t < extremeT) { 151 // we are at the backward end of the history 152 153 // check if a new rough root has been crossed 154 final Transformer previous = transformers[0]; 155 final Transformer next = filter.selectTransformer(previous, rawG, forward); 156 if (next != previous) { 157 // there is a root somewhere between extremeT and t. 158 // the new transformer is valid for t (this is how we have just computed 159 // it above), but it is in fact valid on both sides of the root, so 160 // it was already valid before t and even up to previous time. We store 161 // the switch at extremeT for safety, to ensure the previous transformer 162 // is not applied too close of the root 163 System.arraycopy(updates, 0, updates, 1, updates.length - 1); 164 System.arraycopy(transformers, 0, transformers, 1, transformers.length - 1); 165 updates[0] = extremeT; 166 transformers[0] = next; 167 } 168 169 extremeT = t; 170 171 // apply the transform 172 return next.transformed(rawG); 173 174 } else { 175 // we are in the middle of the history 176 177 // select the transformer 178 for (int i = 0; i < updates.length - 1; ++i) { 179 if (t <= updates[i]) { 180 // apply the transform 181 return transformers[i].transformed(rawG); 182 } 183 } 184 185 return transformers[updates.length - 1].transformed(rawG); 186 187 } 188 } 189 190 } 191 192 /** {@inheritDoc} */ 193 public Action eventOccurred(double t, double[] y, boolean increasing) { 194 // delegate to raw handler, fixing increasing status on the fly 195 return rawHandler.eventOccurred(t, y, filter.getTriggeredIncreasing()); 196 } 197 198 /** {@inheritDoc} */ 199 public void resetState(double t, double[] y) { 200 // delegate to raw handler 201 rawHandler.resetState(t, y); 202 } 203 204}