1 /*
2 * Licensed to the Hipparchus project under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The Hipparchus project licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * https://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.hipparchus.analysis.polynomials;
18
19 import org.hipparchus.CalculusFieldElement;
20 import org.hipparchus.Field;
21 import org.hipparchus.exception.LocalizedCoreFormats;
22 import org.hipparchus.exception.MathIllegalArgumentException;
23 import org.hipparchus.exception.NullArgumentException;
24 import org.hipparchus.util.MathArrays;
25
26 /**
27 * Smoothstep function factory.
28 * <p>
29 * It allows for quick creation of common and generic smoothstep functions as defined
30 * <a href="https://en.wikipedia.org/wiki/Smoothstep">here</a>.
31 */
32 public class SmoothStepFactory {
33
34 /**
35 * Private constructor.
36 * <p>
37 * This class is a utility class, it should neither have a public nor a default constructor. This private constructor
38 * prevents the compiler from generating one automatically.
39 */
40 private SmoothStepFactory() {
41 // Empty constructor
42 }
43
44 /**
45 * Get the {@link SmoothStepFunction clamping smoothstep function}.
46 *
47 * @return clamping smoothstep function
48 */
49 public static SmoothStepFunction getClamp() {
50 return getGeneralOrder(0);
51 }
52
53 /**
54 * Get the {@link SmoothStepFunction quadratic smoothstep function}.
55 *
56 * @return clamping smoothstep function
57 */
58 public static SmoothStepFunction getQuadratic() {
59 // Use a default double array as it will not matter anyway
60 return new QuadraticSmoothStepFunction(new double[] { 0 });
61 }
62
63 /**
64 * Get the {@link SmoothStepFunction cubic smoothstep function}.
65 *
66 * @return cubic smoothstep function
67 */
68 public static SmoothStepFunction getCubic() {
69 return getGeneralOrder(1);
70 }
71
72 /**
73 * Get the {@link SmoothStepFunction quintic smoothstep function}.
74 *
75 * @return quintic smoothstep function
76 */
77 public static SmoothStepFunction getQuintic() {
78 return getGeneralOrder(2);
79 }
80
81 /**
82 * Get the {@link SmoothStepFunction clamping smoothstep function}.
83 *
84 * @param <T> type of the field element
85 * @param field field of the element
86 *
87 * @return clamping smoothstep function
88 */
89 public static <T extends CalculusFieldElement<T>> FieldSmoothStepFunction<T> getClamp(final Field<T> field) {
90 return getFieldGeneralOrder(field, 0);
91 }
92
93 /**
94 * Get the {@link SmoothStepFunction quadratic smoothstep function}.
95 *
96 * @param <T> type of the field element
97 * @param field field of the element
98 *
99 * @return clamping smoothstep function
100 */
101 public static <T extends CalculusFieldElement<T>> FieldSmoothStepFunction<T> getQuadratic(final Field<T> field) {
102 final T[] tempArray = MathArrays.buildArray(field, 1);
103 return new FieldQuadraticSmoothStepFunction<>(tempArray);
104 }
105
106 /**
107 * Get the {@link SmoothStepFunction cubic smoothstep function}.
108 *
109 * @param <T> type of the field element
110 * @param field field of the element
111 *
112 * @return cubic smoothstep function
113 */
114 public static <T extends CalculusFieldElement<T>> FieldSmoothStepFunction<T> getCubic(final Field<T> field) {
115 return getFieldGeneralOrder(field, 1);
116 }
117
118 /**
119 * Get the {@link SmoothStepFunction quintic smoothstep function}.
120 *
121 * @param <T> type of the field element
122 * @param field field of the element
123 *
124 * @return quintic smoothstep function
125 */
126 public static <T extends CalculusFieldElement<T>> FieldSmoothStepFunction<T> getQuintic(final Field<T> field) {
127 return getFieldGeneralOrder(field, 2);
128 }
129
130 /**
131 * Create a {@link SmoothStepFunction smoothstep function} of order <b>2N + 1</b>.
132 * <p>
133 * It uses the general smoothstep equation presented <a href="https://en.wikipedia.org/wiki/Smoothstep">here</a> :
134 * $S_{N}(x) = \sum_{n=0}^{N} \begin{pmatrix} -N-1 \\ n \end{pmatrix} \begin{pmatrix} 2N+1 \\ N-n \end{pmatrix}
135 * x^{N+n+1}$
136 *
137 * @param N determines the order of the output smoothstep function (=2N + 1)
138 *
139 * @return smoothstep function of order <b>2N + 1</b>
140 */
141 public static SmoothStepFunction getGeneralOrder(final int N) {
142
143 final int twoNPlusOne = 2 * N + 1;
144
145 final double[] coefficients = new double[twoNPlusOne + 1];
146
147 int n = N;
148 for (int i = twoNPlusOne; i > N; i--) {
149 coefficients[i] = pascalTriangle(-N - 1, n) * pascalTriangle(2 * N + 1, N - n);
150 n--;
151 }
152
153 return new SmoothStepFunction(coefficients);
154 }
155
156 /**
157 * Create a {@link SmoothStepFunction smoothstep function} of order <b>2N + 1</b>.
158 * <p>
159 * It uses the general smoothstep equation presented <a href="https://en.wikipedia.org/wiki/Smoothstep">here</a> :
160 * $S_{N}(x) = \sum_{n=0}^{N} \begin{pmatrix} -N-1 \\ n \end{pmatrix} \begin{pmatrix} 2N+1 \\ N-n \end{pmatrix}
161 * x^{N+n+1}$
162 *
163 * @param <T> type of the field element
164 * @param field field of the element
165 * @param N determines the order of the output smoothstep function (=2N + 1)
166 *
167 * @return smoothstep function of order <b>2N + 1</b>
168 */
169 public static <T extends CalculusFieldElement<T>> FieldSmoothStepFunction<T> getFieldGeneralOrder(final Field<T> field,
170 final int N) {
171
172 final int twoNPlusOne = 2 * N + 1;
173
174 final T[] coefficients = MathArrays.buildArray(field, twoNPlusOne + 1);
175
176 final T one = field.getOne();
177 int n = N;
178 for (int i = twoNPlusOne; i > N; i--) {
179 coefficients[i] = one.newInstance(pascalTriangle(-N - 1, n) * pascalTriangle(2 * N + 1, N - n));
180 n--;
181 }
182
183 return new FieldSmoothStepFunction<>(coefficients);
184 }
185
186 /**
187 * Returns binomial coefficient without explicit use of factorials, which can't be used with negative integers
188 *
189 * @param k subset in set
190 * @param n set
191 *
192 * @return number of subset {@code k} in global set {@code n}
193 */
194 private static int pascalTriangle(final int k, final int n) {
195
196 int result = 1;
197 for (int i = 0; i < n; i++) {
198 result *= (k - i) / (i + 1);
199 }
200
201 return result;
202 }
203
204 /**
205 * Check that input is between [0:1].
206 *
207 * @param input input to be checked
208 *
209 * @throws MathIllegalArgumentException if input is not between [0:1]
210 */
211 public static void checkBetweenZeroAndOneIncluded(final double input) throws MathIllegalArgumentException {
212 if (input < 0 || input > 1) {
213 throw new MathIllegalArgumentException(
214 LocalizedCoreFormats.INPUT_EXPECTED_BETWEEN_ZERO_AND_ONE_INCLUDED);
215 }
216 }
217
218 /**
219 * Smoothstep function as defined <a href="https://en.wikipedia.org/wiki/Smoothstep">here</a>.
220 * <p>
221 * It is used to do a smooth transition between the "left edge" and the "right edge" with left edge assumed to be smaller
222 * than right edge.
223 * <p>
224 * By definition, for order n greater than 1 and input x, a smoothstep function respects at least the following properties :
225 * <ul>
226 * <li>f(x <= leftEdge) = 0 and f(x >= rightEdge) = 1</li>
227 * <li>f'(leftEdge) = f'(rightEdge) = 0</li>
228 * </ul>
229 * If x is normalized between edges, we have at least :
230 * <ul>
231 * <li>f(x <= 0) = 0 and f(x >= 1) = 1</li>
232 * <li>f'(0) = f'(1) = 0</li>
233 * </ul>
234 * Smoothstep functions of higher order n will have their higher time derivatives also equal to zero at edges...
235 */
236 public static class SmoothStepFunction extends PolynomialFunction {
237
238 /** Serializable UID. */
239 private static final long serialVersionUID = 20230113L;
240
241 /**
242 * Construct a smoothstep with the given coefficients. The first element of the coefficients array is the constant
243 * term. Higher degree coefficients follow in sequence. The degree of the resulting polynomial is the index of the
244 * last non-null element of the array, or 0 if all elements are null.
245 * <p>
246 * The constructor makes a copy of the input array and assigns the copy to the coefficients property.</p>
247 *
248 * @param c Smoothstep polynomial coefficients.
249 *
250 * @throws NullArgumentException if {@code c} is {@code null}.
251 * @throws MathIllegalArgumentException if {@code c} is empty.
252 */
253 private SmoothStepFunction(final double[] c) throws MathIllegalArgumentException, NullArgumentException {
254 super(c);
255 }
256
257 /**
258 * Compute the value of the smoothstep for the given argument normalized between edges.
259 *
260 * @param xNormalized Normalized argument for which the function value should be computed. It is expected to be
261 * between [0:1] and will throw an exception otherwise.
262 *
263 * @return the value of the polynomial at the given point.
264 *
265 * @see org.hipparchus.analysis.UnivariateFunction#value(double)
266 */
267 @Override
268 public double value(final double xNormalized) {
269 checkBetweenZeroAndOneIncluded(xNormalized);
270 return super.value(xNormalized);
271 }
272
273 /**
274 * Compute the value of the smoothstep function for the given edges and argument.
275 * <p>
276 * Note that right edge is expected to be greater than left edge. It will throw an exception otherwise.
277 *
278 * @param leftEdge left edge
279 * @param rightEdge right edge
280 * @param x Argument for which the function value should be computed
281 *
282 * @return the value of the polynomial at the given point
283 *
284 * @throws MathIllegalArgumentException if right edge is greater than left edge
285 * @see org.hipparchus.analysis.UnivariateFunction#value(double)
286 */
287 public double value(final double leftEdge, final double rightEdge, final double x)
288 throws MathIllegalArgumentException {
289
290 checkInputEdges(leftEdge, rightEdge);
291
292 final double xClamped = clampInput(leftEdge, rightEdge, x);
293
294 final double xNormalized = normalizeInput(leftEdge, rightEdge, xClamped);
295
296 return super.value(xNormalized);
297 }
298
299 /**
300 * Check that left edge is lower than right edge. Otherwise, throw an exception.
301 *
302 * @param leftEdge left edge
303 * @param rightEdge right edge
304 */
305 protected void checkInputEdges(final double leftEdge, final double rightEdge) {
306 if (leftEdge > rightEdge) {
307 throw new MathIllegalArgumentException(LocalizedCoreFormats.RIGHT_EDGE_GREATER_THAN_LEFT_EDGE,
308 leftEdge, rightEdge);
309 }
310 }
311
312 /**
313 * Clamp input between edges.
314 *
315 * @param leftEdge left edge
316 * @param rightEdge right edge
317 * @param x input to clamp
318 *
319 * @return clamped input
320 */
321 protected double clampInput(final double leftEdge, final double rightEdge, final double x) {
322 if (x <= leftEdge) {
323 return leftEdge;
324 }
325 if (x >= rightEdge) {
326 return rightEdge;
327 }
328 return x;
329 }
330
331 /**
332 * Normalize input between left and right edges.
333 *
334 * @param leftEdge left edge
335 * @param rightEdge right edge
336 * @param x input to normalize
337 *
338 * @return normalized input
339 */
340 protected double normalizeInput(final double leftEdge, final double rightEdge, final double x) {
341 return (x - leftEdge) / (rightEdge - leftEdge);
342 }
343 }
344
345 /**
346 * Specific smoothstep function that cannot be built using the {@link #getGeneralOrder(int)}.
347 * <p>
348 * Methods inherited from {@link PolynomialFunction} <em>should not be used</em> as they will not be true to the actual
349 * function.
350 *
351 * @see PolynomialFunction
352 */
353 public static class QuadraticSmoothStepFunction extends SmoothStepFunction {
354
355 /** Serializable UID. */
356 private static final long serialVersionUID = 20230422L;
357
358 /**
359 * Construct a smoothstep with the given coefficients. The first element of the coefficients array is the constant
360 * term. Higher degree coefficients follow in sequence. The degree of the resulting polynomial is the index of the
361 * last non-null element of the array, or 0 if all elements are null.
362 * <p>
363 * The constructor makes a copy of the input array and assigns the copy to the coefficients property.</p>
364 *
365 * @param c Smoothstep polynomial coefficients.
366 *
367 * @throws NullArgumentException if {@code c} is {@code null}.
368 * @throws MathIllegalArgumentException if {@code c} is empty.
369 */
370 private QuadraticSmoothStepFunction(final double[] c) throws MathIllegalArgumentException, NullArgumentException {
371 super(c);
372 }
373
374 /**
375 * Compute the value of the smoothstep function for the given edges and argument.
376 * <p>
377 * Note that right edge is expected to be greater than left edge. It will throw an exception otherwise.
378 *
379 * @param leftEdge left edge
380 * @param rightEdge right edge
381 * @param x Argument for which the function value should be computed
382 *
383 * @return the value of the polynomial at the given point
384 *
385 * @throws MathIllegalArgumentException if right edge is greater than left edge
386 * @see org.hipparchus.analysis.UnivariateFunction#value(double)
387 */
388 @Override
389 public double value(final double leftEdge, final double rightEdge, final double x)
390 throws MathIllegalArgumentException {
391
392 checkInputEdges(leftEdge, rightEdge);
393
394 final double xClamped = clampInput(leftEdge, rightEdge, x);
395
396 final double xNormalized = normalizeInput(leftEdge, rightEdge, xClamped);
397
398 return value(xNormalized);
399 }
400
401 /**
402 * Compute the value of the quadratic smoothstep for the given argument normalized between edges.
403 *
404 * @param xNormalized Normalized argument for which the function value should be computed. It is expected to be
405 * between [0:1] and will throw an exception otherwise.
406 *
407 * @return the value of the polynomial at the given point.
408 *
409 * @see org.hipparchus.analysis.UnivariateFunction#value(double)
410 */
411 @Override
412 public double value(final double xNormalized) {
413 checkBetweenZeroAndOneIncluded(xNormalized);
414
415 if (xNormalized >= 0 && xNormalized <= 0.5) {
416 return 2 * xNormalized * xNormalized;
417 }
418 else {
419 return 4 * xNormalized - 2 * xNormalized * xNormalized - 1;
420 }
421 }
422 }
423
424 /**
425 * Smoothstep function as defined <a href="https://en.wikipedia.org/wiki/Smoothstep">here</a>.
426 * <p>
427 * It is used to do a smooth transition between the "left edge" and the "right edge" with left edge assumed to be smaller
428 * than right edge.
429 * <p>
430 * By definition, for order n greater than 1 and input x, a smoothstep function respects at least the following properties :
431 * <ul>
432 * <li>f(x <= leftEdge) = 0 and f(x >= rightEdge) = 1</li>
433 * <li>f'(leftEdge) = f'(rightEdge) = 0</li>
434 * </ul>
435 * If x is normalized between edges, we have at least :
436 * <ul>
437 * <li>f(x <= 0) = 0 and f(x >= 1) = 1</li>
438 * <li>f'(0) = f'(1) = 0</li>
439 * </ul>
440 * Smoothstep functions of higher order n will have their higher time derivatives also equal to zero at edges...
441 *
442 * @param <T> type of the field element
443 */
444 public static class FieldSmoothStepFunction<T extends CalculusFieldElement<T>> extends FieldPolynomialFunction<T> {
445
446 /**
447 * Construct a smoothstep with the given coefficients. The first element of the coefficients array is the constant
448 * term. Higher degree coefficients follow in sequence. The degree of the resulting polynomial is the index of the
449 * last non-null element of the array, or 0 if all elements are null.
450 * <p>
451 * The constructor makes a copy of the input array and assigns the copy to the coefficients property.</p>
452 *
453 * @param c Smoothstep polynomial coefficients.
454 *
455 * @throws NullArgumentException if {@code c} is {@code null}.
456 * @throws MathIllegalArgumentException if {@code c} is empty.
457 */
458 private FieldSmoothStepFunction(final T[] c) throws MathIllegalArgumentException, NullArgumentException {
459 super(c);
460 }
461
462 /**
463 * Compute the value of the smoothstep for the given argument normalized between edges.
464 *
465 * @param xNormalized Normalized argument for which the function value should be computed. It is expected to be
466 * between [0:1] and will throw an exception otherwise.
467 *
468 * @return the value of the polynomial at the given point.
469 *
470 * @see org.hipparchus.analysis.UnivariateFunction#value(double)
471 */
472 @Override
473 public T value(final double xNormalized) {
474 checkBetweenZeroAndOneIncluded(xNormalized);
475 return super.value(xNormalized);
476 }
477
478 /**
479 * Compute the value of the smoothstep for the given argument normalized between edges.
480 *
481 * @param xNormalized Normalized argument for which the function value should be computed. It is expected to be
482 * between [0:1] and will throw an exception otherwise.
483 *
484 * @return the value of the polynomial at the given point.
485 *
486 * @see org.hipparchus.analysis.UnivariateFunction#value(double)
487 */
488 @Override
489 public T value(final T xNormalized) {
490 checkBetweenZeroAndOneIncluded(xNormalized.getReal());
491 return super.value(xNormalized);
492 }
493
494 /**
495 * Compute the value of the smoothstep function for the given edges and argument.
496 * <p>
497 * Note that right edge is expected to be greater than left edge. It will throw an exception otherwise.
498 *
499 * @param leftEdge left edge
500 * @param rightEdge right edge
501 * @param x Argument for which the function value should be computed
502 *
503 * @return the value of the polynomial at the given point
504 *
505 * @throws MathIllegalArgumentException if right edge is greater than left edge
506 * @see org.hipparchus.analysis.UnivariateFunction#value(double)
507 */
508 public T value(final double leftEdge, final double rightEdge, final T x)
509 throws MathIllegalArgumentException {
510
511 checkInputEdges(leftEdge, rightEdge);
512
513 final T xClamped = clampInput(leftEdge, rightEdge, x);
514
515 final T xNormalized = normalizeInput(leftEdge, rightEdge, xClamped);
516
517 return super.value(xNormalized);
518 }
519
520 /**
521 * Check that left edge is lower than right edge. Otherwise, throw an exception.
522 *
523 * @param leftEdge left edge
524 * @param rightEdge right edge
525 */
526 protected void checkInputEdges(final double leftEdge, final double rightEdge) {
527 if (leftEdge > rightEdge) {
528 throw new MathIllegalArgumentException(LocalizedCoreFormats.RIGHT_EDGE_GREATER_THAN_LEFT_EDGE,
529 leftEdge, rightEdge);
530 }
531 }
532
533 /**
534 * Clamp input between edges.
535 *
536 * @param leftEdge left edge
537 * @param rightEdge right edge
538 * @param x input to clamp
539 *
540 * @return clamped input
541 */
542 protected T clampInput(final double leftEdge, final double rightEdge, final T x) {
543 if (x.getReal() <= leftEdge) {
544 return x.getField().getOne().newInstance(leftEdge);
545 }
546 if (x.getReal() >= rightEdge) {
547 return x.getField().getOne().newInstance(rightEdge);
548 }
549 return x;
550 }
551
552 /**
553 * Normalize input between left and right edges.
554 *
555 * @param leftEdge left edge
556 * @param rightEdge right edge
557 * @param x input to normalize
558 *
559 * @return normalized input
560 */
561 protected T normalizeInput(final double leftEdge, final double rightEdge, final T x) {
562 return x.subtract(leftEdge).divide(rightEdge - leftEdge);
563 }
564 }
565
566 /**
567 * Specific smoothstep function that cannot be built using the {@link #getGeneralOrder(int)}.
568 * <p>
569 * Methods inherited from {@link PolynomialFunction} <em>should not be used</em> as they will not be true to the actual
570 * function.
571 *
572 * @param <T> type of the field element
573 *
574 * @see PolynomialFunction
575 */
576 private static class FieldQuadraticSmoothStepFunction<T extends CalculusFieldElement<T>>
577 extends FieldSmoothStepFunction<T> {
578
579 /**
580 * Construct a smoothstep with the given coefficients. The first element of the coefficients array is the constant
581 * term. Higher degree coefficients follow in sequence. The degree of the resulting polynomial is the index of the
582 * last non-null element of the array, or 0 if all elements are null.
583 * <p>
584 * The constructor makes a copy of the input array and assigns the copy to the coefficients property.</p>
585 *
586 * @param c Smoothstep polynomial coefficients.
587 *
588 * @throws NullArgumentException if {@code c} is {@code null}.
589 * @throws MathIllegalArgumentException if {@code c} is empty.
590 */
591 private FieldQuadraticSmoothStepFunction(final T[] c) throws MathIllegalArgumentException, NullArgumentException {
592 super(c);
593 }
594
595 /**
596 * Compute the value of the smoothstep function for the given edges and argument.
597 * <p>
598 * Note that right edge is expected to be greater than left edge. It will throw an exception otherwise.
599 *
600 * @param leftEdge left edge
601 * @param rightEdge right edge
602 * @param x Argument for which the function value should be computed
603 *
604 * @return the value of the polynomial at the given point
605 *
606 * @throws MathIllegalArgumentException if right edge is greater than left edge
607 * @see org.hipparchus.analysis.UnivariateFunction#value(double)
608 */
609 @Override
610 public T value(final double leftEdge, final double rightEdge, final T x)
611 throws MathIllegalArgumentException {
612
613 checkInputEdges(leftEdge, rightEdge);
614
615 final T xClamped = clampInput(leftEdge, rightEdge, x);
616
617 final T xNormalized = normalizeInput(leftEdge, rightEdge, xClamped);
618
619 return value(xNormalized);
620 }
621
622 /**
623 * Compute the value of the quadratic smoothstep for the given argument normalized between edges.
624 * <p>
625 *
626 * @param xNormalized Normalized argument for which the function value should be computed. It is expected to be
627 * between [0:1] and will throw an exception otherwise.
628 *
629 * @return the value of the polynomial at the given point.
630 *
631 * @see org.hipparchus.analysis.UnivariateFunction#value(double)
632 */
633 @Override
634 public T value(final double xNormalized) {
635 checkBetweenZeroAndOneIncluded(xNormalized);
636
637 final Field<T> field = getField();
638 final T one = field.getOne();
639
640 if (xNormalized >= 0 && xNormalized <= 0.5) {
641 return one.newInstance(2. * xNormalized * xNormalized);
642 }
643 else {
644 return one.newInstance(4. * xNormalized - 2. * xNormalized * xNormalized - 1.);
645 }
646 }
647
648 /**
649 * Compute the value of the quadratic smoothstep for the given argument normalized between edges.
650 * <p>
651 *
652 * @param xNormalized Normalized argument for which the function value should be computed. It is expected to be
653 * between [0:1] and will throw an exception otherwise.
654 *
655 * @return the value of the polynomial at the given point.
656 *
657 * @see org.hipparchus.analysis.UnivariateFunction#value(double)
658 */
659 @Override
660 public T value(final T xNormalized) {
661 checkBetweenZeroAndOneIncluded(xNormalized.getReal());
662
663 if (xNormalized.getReal() >= 0 && xNormalized.getReal() <= 0.5) {
664 return xNormalized.multiply(xNormalized).multiply(2.);
665 }
666 else {
667 final T one = getField().getOne();
668 return one.linearCombination(4., xNormalized, -2., xNormalized.multiply(xNormalized)).subtract(1.);
669 }
670 }
671 }
672
673 }