Repository: sdpython/onnxcustom (language: python)

Issue:
Implements classificiation loss



Discussion / hints:

Fix patch:
diff --git a/.circleci/config.yml b/.circleci/config.yml
index 52302045..3488e6c2 100644
--- a/.circleci/config.yml
+++ b/.circleci/config.yml
@@ -49,8 +49,8 @@ jobs:
           name: flake8
           command: |
             . venv/bin/activate
-            python -m flake8 onnxcustom --max-line-length 90
-            python -m flake8 _doc/examples --max-line-length 90
+            python -m flake8 onnxcustom --max-line-length 100
+            python -m flake8 _doc/examples --max-line-length 100
 
       - run:
           name: run tests
diff --git a/.gitignore b/.gitignore
index 0ceeffbd..1cc45db6 100644
--- a/.gitignore
+++ b/.gitignore
@@ -60,3 +60,4 @@ _unittests/ut_documentation/plot_linear_regression.png
 _unittests/ut_documentation/summary.csv
 _unittests/ut_documentation/_test_example.txt
 _unittests/ut_documentation/_test_example.txt
+something
diff --git a/_doc/examples/plot_funny_sigmoid.py b/_doc/examples/plot_funny_sigmoid.py
new file mode 100644
index 00000000..f220ab05
--- /dev/null
+++ b/_doc/examples/plot_funny_sigmoid.py
@@ -0,0 +1,160 @@
+"""
+.. _l-example-discrepencies-sigmoid:
+
+Funny discrepancies
+===================
+
+Function sigmoid is :math:`sig(x) = \\frac{1}{1 + e^{-x}}`.
+For small or high value, implementation has to do approximation
+and they are not always the same. It may be a tradeoff between
+precision and computation time...
+It is always a tradeoff.
+
+.. index:: discrepencies, sigmoid
+
+.. contents::
+    :local:
+
+
+Precision
++++++++++
+
+This section compares the precision of a couple of implementations
+of the ssigmoid function. The custom implementation is done with
+a Taylor expansion of exponential function:
+:math:`e^x \\sim 1 + x + \\frac{x^2}{2} + ... + \\frac{x^n}{n!} + o(x^n)`.
+
+"""
+import time
+import numpy
+import pandas
+from tqdm import tqdm
+from scipy.special import expit
+
+from skl2onnx.algebra.onnx_ops import OnnxSigmoid
+from skl2onnx.common.data_types import FloatTensorType
+from onnxruntime import InferenceSession
+from mlprodict.onnxrt import OnnxInference
+from onnxcustom import get_max_opset
+import matplotlib.pyplot as plt
+
+one = numpy.array([1], dtype=numpy.float64)
+
+
+def taylor_approximation_exp(x, degre=50):
+    y = numpy.zeros(x.shape, dtype=x.dtype)
+    a = numpy.ones(x.shape, dtype=x.dtype)
+    for i in range(1, degre + 1):
+        a *= x / i
+        y += a
+    return y
+
+
+def taylor_sigmoid(x, degre=50):
+    den = one + taylor_approximation_exp(-x, degre)
+    return one / (den)
+
+
+opset = get_max_opset()
+N = 300
+min_values = [-20 + float(i) * 10 / N for i in range(N)]
+data = numpy.array([0], dtype=numpy.float32)
+
+node = OnnxSigmoid('X', op_version=opset, output_names=['Y'])
+onx = node.to_onnx({'X': FloatTensorType()},
+                   {'Y': FloatTensorType()},
+                   target_opset=opset)
+rts = ['numpy', 'python', 'onnxruntime', 'taylor20', 'taylor40']
+
+oinf = OnnxInference(onx)
+sess = InferenceSession(onx.SerializeToString())
+
+graph = []
+for mv in tqdm(min_values):
+    data[0] = mv
+    for rt in rts:
+        lab = ""
+        if rt == 'numpy':
+            y = expit(data)
+        elif rt == 'python':
+            y = oinf.run({'X': data})['Y']
+            # * 1.2 to avoid curves to be superimposed
+            y *= 1.2
+            lab = "x1.2"
+        elif rt == 'onnxruntime':
+            y = sess.run(None, {'X': data})[0]
+        elif rt == 'taylor40':
+            y = taylor_sigmoid(data, 40)
+            # * 0.8 to avoid curves to be superimposed
+            y *= 0.8
+            lab = "x0.8"
+        elif rt == 'taylor20':
+            y = taylor_sigmoid(data, 20)
+            # * 0.6 to avoid curves to be superimposed
+            y *= 0.6
+            lab = "x0.6"
+        else:
+            raise AssertionError("Unknown runtime %r." % rt)
+        value = y[0]
+        graph.append(dict(rt=rt + lab, x=mv, y=value))
+
+#############################################
+# Graph.
+
+_, ax = plt.subplots(1, 1, figsize=(12, 4))
+df = pandas.DataFrame(graph)
+piv = df.pivot('x', 'rt', 'y')
+print(piv.T.head())
+piv.plot(ax=ax, logy=True)
+
+##############################################
+# :math:`log(sig(x)) = -log(1 + e^{-x})`. When *x* is very negative,
+# :math:`log(sig(x)) \\sim -x`. That explains the graph.
+# We also see :epkg:`onnxruntime` is less precise for these values.
+# What's the benefit?
+#
+# Computation time
+# ++++++++++++++++
+
+graph = []
+for mv in tqdm(min_values):
+    data = numpy.array([mv] * 10000, dtype=numpy.float32)
+    for rt in rts:
+        begin = time.perf_counter()
+        if rt == 'numpy':
+            y = expit(data)
+        elif rt == 'python':
+            y = oinf.run({'X': data})['Y']
+        elif rt == 'onnxruntime':
+            y = sess.run(None, {'X': data})[0]
+        elif rt == 'taylor40':
+            y = taylor_sigmoid(data, 40)
+        elif rt == 'taylor20':
+            y = taylor_sigmoid(data, 20)
+        else:
+            raise AssertionError("Unknown runtime %r." % rt)
+        duration = time.perf_counter() - begin
+        graph.append(dict(rt=rt, x=mv, y=duration))
+
+_, ax = plt.subplots(1, 1, figsize=(12, 4))
+df = pandas.DataFrame(graph)
+piv = df.pivot('x', 'rt', 'y')
+piv.plot(ax=ax, logy=True)
+
+#############################################
+# Conclusion
+# ++++++++++
+#
+# The implementation from :epkg:`onnxruntime` is faster but
+# is much less contiguous for extremes. It explains why
+# probabilities may be much different when an observation
+# is far from every classification border. In that case,
+# :epkg:`onnxruntime` implementation of the sigmoid function
+# returns 0 when :func:`numpy.sigmoid` returns a smoother value.
+# Probabilites of logistic regression are obtained after the raw
+# scores are transformed with the sigmoid function and
+# normalized. If the raw scores are very negative,
+# the sum of probabilities becomes null with :epkg:`onnxruntime`.
+# The normalization fails.
+
+# plt.show()
diff --git a/_doc/sphinxdoc/source/api/training.rst b/_doc/sphinxdoc/source/api/training.rst
index 63eae5b4..8bd39d19 100644
--- a/_doc/sphinxdoc/source/api/training.rst
+++ b/_doc/sphinxdoc/source/api/training.rst
@@ -26,7 +26,7 @@ BaseEstimator
 
 Ancestor to both classes wrapping :epkg:`onnxruntime` API.
 
-.. autosignature:: onnxcustom.training.base_estimator.BaseEstimator
+.. autosignature:: onnxcustom.training._base_estimator.BaseEstimator
     :members:
 
 Exceptions
@@ -52,6 +52,8 @@ be combination of L1, L2 losses and L1, L2 penalties.
 
 .. autosignature:: onnxcustom.utils.orttraining_helper.get_train_initializer
 
+.. autosignature:: onnxcustom.utils.onnx_rewriter.onnx_rewrite_operator
+
 .. _l-api-prt-gradient-optimizer:
 
 OrtGradientOptimizer
@@ -100,11 +102,14 @@ LearningLoss
 .. autosignature:: onnxcustom.training.sgd_learning_loss.ElasticLearningLoss
     :members:
 
+.. autosignature:: onnxcustom.training.sgd_learning_loss.NegLogLearningLoss
+    :members:
+
 .. autosignature:: onnxcustom.training.sgd_learning_loss.SquareLearningLoss
     :members:
 
-Loss function
-+++++++++++++
+Loss functions
+++++++++++++++
 
 .. autosignature:: onnxcustom.utils.onnx_function.function_onnx_graph
 
diff --git a/_doc/sphinxdoc/source/tutorial_onnxruntime/training_ort_api.rst b/_doc/sphinxdoc/source/tutorial_onnxruntime/training_ort_api.rst
index a2737bc3..2b1440eb 100644
--- a/_doc/sphinxdoc/source/tutorial_onnxruntime/training_ort_api.rst
+++ b/_doc/sphinxdoc/source/tutorial_onnxruntime/training_ort_api.rst
@@ -53,7 +53,7 @@ optimized weights.
 
 :epkg:`onnxruntime-training` does not implement loss functions.
 That must be done independently. That's what function
-:func:`onnxcustom.utils.orttraining_helper.add_loss_output>` does.
+:func:`onnxcustom.utils.orttraining_helper.add_loss_output` does.
 It implements a couple of usual losses in ONNX.
 Another function :func:`onnxcustom.utils.orttraining_helper.get_train_initializer`
 guesses all the coefficients of an ONNX graph if the user does not specify any.
diff --git a/_doc/sphinxdoc/source/tutorial_skl/tutorial_1_simple.rst b/_doc/sphinxdoc/source/tutorial_skl/tutorial_1_simple.rst
index 24b7df13..be732f73 100644
--- a/_doc/sphinxdoc/source/tutorial_skl/tutorial_1_simple.rst
+++ b/_doc/sphinxdoc/source/tutorial_skl/tutorial_1_simple.rst
@@ -21,6 +21,7 @@ used in the ONNX graph.
     ../gyexamples/plot_dbegin_options_list
     ../gyexamples/plot_dbegin_options_zipmap
     ../gyexamples/plot_ebegin_float_double
+    ../gyexamples/plot_funny_sigmoid
     ../gyexamples/plot_fbegin_investigate
     ../gyexamples/plot_gbegin_dataframe
     ../gyexamples/plot_gbegin_transfer_learning
diff --git a/_doc/sphinxdoc/source/tutorial_training/tutorial_6_training_partial.rst b/_doc/sphinxdoc/source/tutorial_training/tutorial_6_training_partial.rst
index 9deb14db..9f24f9a8 100644
--- a/_doc/sphinxdoc/source/tutorial_training/tutorial_6_training_partial.rst
+++ b/_doc/sphinxdoc/source/tutorial_training/tutorial_6_training_partial.rst
@@ -117,11 +117,11 @@ Cache
 +++++
 
 Base class :class:`BaseLearningOnnx
-<onnxcustom.training.base_onnx_function.BaseLearningOnnx>` implements
+<onnxcustom.training._base_onnx_function.BaseLearningOnnx>` implements
 methods :meth:`_bind_input_ortvalue
-<onnxcustom.training.base_onnx_function.BaseLearningOnnx._bind_input_ortvalue>`
+<onnxcustom.training._base_onnx_function.BaseLearningOnnx._bind_input_ortvalue>`
 and :meth:`_bind_output_ortvalue
-<onnxcustom.training.base_onnx_function.BaseLearningOnnx._bind_output_ortvalue>`
+<onnxcustom.training._base_onnx_function.BaseLearningOnnx._bind_output_ortvalue>`
 used by the three components mentioned above. They cache the binded pointers
 (the value returns by `c_ortvalue.data_ptr()` and do not bind again
 if the method is called again with a different `OrtValue` but a same pointer
diff --git a/appveyor.yml b/appveyor.yml
index 674e9454..474def38 100644
--- a/appveyor.yml
+++ b/appveyor.yml
@@ -11,9 +11,9 @@ init:
 install:
   - "%PYTHON%\\python -m pip install -r requirements.txt"
   - "%PYTHON%\\python -m pip install -r requirements-dev.txt"
-  - "%PYTHON%\\python -m pip install http://www.xavierdupre.fr/enseignement/setup/onnxruntime_training-1.10.992+cpu-cp39-cp39-win_amd64.whl"
-  - "%PYTHON%\\python -m pip install http://www.xavierdupre.fr/enseignement/setup/torch_interop_utils-1.10.0+cpu-cp39-cp39-win_amd64.whl"
-  - "%PYTHON%\\python -m pip install http://www.xavierdupre.fr/enseignement/setup/aten_op_executor-1.10.0+cpu-cp39-cp39-win_amd64.whl"
+  - "%PYTHON%\\python -m pip install http://www.xavierdupre.fr/enseignement/setup/onnxruntime_training-1.11.993+cpu-cp39-cp39-win_amd64.whl"
+  - "%PYTHON%\\python -m pip install http://www.xavierdupre.fr/enseignement/setup/torch_interop_utils-1.10.1+cpu-cp39-cp39-win_amd64.whl"
+  - "%PYTHON%\\python -m pip install http://www.xavierdupre.fr/enseignement/setup/aten_op_executor-1.10.1+cpu-cp39-cp39-win_amd64.whl"
 
 build: off
 
@@ -22,8 +22,8 @@ before_test:
 
 test_script:
   - "%PYTHON%\\python setup.py unittests"
-  - "%PYTHON%\\python -m flake8 onnxcustom --max-line-length 90"
-  - "%PYTHON%\\python -m flake8  _doc/examples --max-line-length 90"
+  - "%PYTHON%\\python -m flake8 onnxcustom --max-line-length 100"
+  - "%PYTHON%\\python -m flake8 _doc/examples --max-line-length 100"
 
 after_test:
   - "%PYTHON%\\python setup.py bdist_wheel"
diff --git a/onnxcustom/training/_base.py b/onnxcustom/training/_base.py
index a73d1b53..ec47c171 100644
--- a/onnxcustom/training/_base.py
+++ b/onnxcustom/training/_base.py
@@ -72,7 +72,7 @@ def walk_through(obj, prefix="", only_name=True):
                         print('+', name)
                     else:
                         name = "%s.%s" % (prefix, k) if prefix else k
-                        print('\n++++++', name)
+                        print('\\n++++++', name)
                         print()
                         bf = io.BytesIO(v)
                         onx = onnx.load(bf)
diff --git a/onnxcustom/training/base_estimator.py b/onnxcustom/training/_base_estimator.py
similarity index 98%
rename from onnxcustom/training/base_estimator.py
rename to onnxcustom/training/_base_estimator.py
index bc26a76d..d41e01a2 100644
--- a/onnxcustom/training/base_estimator.py
+++ b/onnxcustom/training/_base_estimator.py
@@ -8,8 +8,8 @@
 from ..utils.onnxruntime_helper import (
     get_ort_device, ort_device_to_string)
 from ._base import BaseOnnxClass
+from ._base_onnx_function import BaseLearningOnnx
 from .sgd_learning_rate import BaseLearningRate
-from .base_onnx_function import BaseLearningOnnx
 
 
 class BaseEstimator(BaseOnnxClass):
diff --git a/onnxcustom/training/base_onnx_function.py b/onnxcustom/training/_base_onnx_function.py
similarity index 100%
rename from onnxcustom/training/base_onnx_function.py
rename to onnxcustom/training/_base_onnx_function.py
diff --git a/onnxcustom/training/optimizers.py b/onnxcustom/training/optimizers.py
index 685964d9..99f64c6d 100644
--- a/onnxcustom/training/optimizers.py
+++ b/onnxcustom/training/optimizers.py
@@ -12,7 +12,7 @@
     numpy_to_ort_value, device_to_providers)
 from .data_loader import OrtDataLoader
 from .excs import ConvergenceError, EvaluationError
-from .base_estimator import BaseEstimator
+from ._base_estimator import BaseEstimator
 
 
 class OrtGradientOptimizer(BaseEstimator):
diff --git a/onnxcustom/training/optimizers_partial.py b/onnxcustom/training/optimizers_partial.py
index eec6b38d..09c8ab8c 100644
--- a/onnxcustom/training/optimizers_partial.py
+++ b/onnxcustom/training/optimizers_partial.py
@@ -14,7 +14,7 @@
 from ..utils.print_helper import str_ortvalue
 from ..utils.orttraining_helper import get_train_initializer
 from .ortgradient import OrtGradientForwardBackward
-from .base_estimator import BaseEstimator
+from ._base_estimator import BaseEstimator
 from .sgd_learning_loss import BaseLearningLoss
 from .sgd_learning_penalty import BaseLearningPenalty
 from .data_loader import OrtDataLoader
diff --git a/onnxcustom/training/ortgradient.py b/onnxcustom/training/ortgradient.py
index 717671b4..9967cc2d 100644
--- a/onnxcustom/training/ortgradient.py
+++ b/onnxcustom/training/ortgradient.py
@@ -3,7 +3,9 @@
 @file
 @brief Gradient with :epkg:`onnxruntime-training` forward backward.
 """
+import os
 import logging
+import warnings
 from io import BytesIO
 import onnx
 from onnx.numpy_helper import to_array
@@ -506,6 +508,40 @@ def __init__(self):
         self.states_ = []
         self.saved_tensors_ = None
 
+    @classmethod
+    def save_onnx_graph(cls, folder, prefix=None, suffix=None):
+        """
+        Saves onnx graph stored in this class.
+        """
+        if folder is None:
+            return None  # pragma: no cover
+        if prefix is None:
+            prefix = ''
+        if suffix is None:
+            suffix = ''
+        if isinstance(folder, str) and not os.path.exists(folder):
+            raise FileNotFoundError(  # pragma: no cover
+                "Folder %r does not exist." % folder)
+        saved = {}
+        for k, v in cls.__dict__.items():
+            if hasattr(v, "SerializeToString"):
+                if isinstance(folder, str):
+                    name = "%s%s%s.%s.onnx" % (
+                        prefix, cls.__name__, suffix, k)
+                    filename = os.path.join(folder, name)
+                    if os.path.exists(filename):
+                        warnings.warn(  # pragma: no cover
+                            "Filename %r already exists." % filename)
+                    with open(filename, "wb") as f:
+                        f.write(v.SerializeToString())
+                    saved[k] = filename
+                else:
+                    saved[k] = v.SerializeToString()
+            elif hasattr(v, "save_onnx_graph"):
+                saved[k] = v.save_onnx_graph(
+                    folder, prefix=prefix, suffix="%s.%s" % (suffix, k))
+        return saved
+
     @staticmethod
     def device_name(device):
         """
diff --git a/onnxcustom/training/sgd_learning_loss.py b/onnxcustom/training/sgd_learning_loss.py
index 2cbdb9eb..a8104f11 100644
--- a/onnxcustom/training/sgd_learning_loss.py
+++ b/onnxcustom/training/sgd_learning_loss.py
@@ -6,7 +6,7 @@
 from onnxruntime import SessionOptions, InferenceSession, RunOptions
 from ..utils.onnx_function import function_onnx_graph
 from ..utils.onnxruntime_helper import device_to_providers
-from .base_onnx_function import BaseLearningOnnx
+from ._base_onnx_function import BaseLearningOnnx
 
 
 class BaseLearningLoss(BaseLearningOnnx):
@@ -76,7 +76,8 @@ def select(class_name, **kwargs):
             return class_name
         cls = {SquareLearningLoss: ['square_error', 'square'],
                AbsoluteLearningLoss: ['absolute_error', 'absolute'],
-               ElasticLearningLoss: ['elastic_error', 'elastic']}
+               ElasticLearningLoss: ['elastic_error', 'elastic'],
+               NegLogLearningLoss: ['log', 'neglog', 'logloss']}
         for cl, aliases in cls.items():
             if class_name == cl.__class__.__name__ or class_name in aliases:
                 return cl(**kwargs)
@@ -171,3 +172,46 @@ def build_onnx_function(self, opset, device, weight_name):
             providers=device_to_providers(device))
         self.loss_grad_sess_bind_ = (
             self.loss_grad_sess_.io_binding()._iobinding)
+
+
+class NegLogLearningLoss(BaseLearningLoss):
+    """
+    Implements a negative log loss
+    `'log(yt, yp) = -(1-yt)\\log(1-yp) - yt\\log(yp)`,
+    this only works for a binary classification where *yp* is the
+    predicted probability, *yt* is the expected probability.
+    *yt* is expected to be binary, *yp* is a matrix with two
+    columns, the sum on every line is 1.
+    However, this loss is usually applied after a function softmax
+    and the gradient is directly computed from the loss to the
+    raw score before they are processed through the softmax function
+    (see class `Log
+    <https://github.com/scikit-learn/scikit-learn/blob/main/sklearn/
+    linear_model/_sgd_fast.pyx#L236>`_).
+
+    :param eps: clipping value for probabilities,
+        avoids computing `log(0)`
+    :param probability_function: function to convert
+        raw scores into probabilities, default value is `sigmoid`
+        for a logistic regression
+    """
+
+    def __init__(self, eps=1e-5, probability_function='sigmoid'):
+        BaseLearningLoss.__init__(self)
+        self.eps = eps
+        self.probability_function = probability_function
+
+    def build_onnx_function(self, opset, device, weight_name):
+        so = SessionOptions()
+        so.log_severity_level = 4
+
+        # loss_grad
+        fct_name = "grad_%s_neg_log_loss_error" % self.probability_function
+        self.loss_grad_onnx_ = function_onnx_graph(
+            fct_name, target_opset=opset,
+            weight_name=weight_name, eps=self.eps)
+        self.loss_grad_sess_ = InferenceSession(
+            self.loss_grad_onnx_.SerializeToString(), so,
+            providers=device_to_providers(device))
+        self.loss_grad_sess_bind_ = (
+            self.loss_grad_sess_.io_binding()._iobinding)
diff --git a/onnxcustom/training/sgd_learning_penalty.py b/onnxcustom/training/sgd_learning_penalty.py
index 3d8bc151..7a8e710e 100644
--- a/onnxcustom/training/sgd_learning_penalty.py
+++ b/onnxcustom/training/sgd_learning_penalty.py
@@ -6,7 +6,7 @@
 from onnxruntime import SessionOptions, InferenceSession, RunOptions
 from ..utils.onnx_function import function_onnx_graph
 from ..utils.onnxruntime_helper import device_to_providers
-from .base_onnx_function import BaseLearningOnnx
+from ._base_onnx_function import BaseLearningOnnx
 
 
 class BaseLearningPenalty(BaseLearningOnnx):
diff --git a/onnxcustom/training/sgd_learning_rate.py b/onnxcustom/training/sgd_learning_rate.py
index ec2db721..8f5185eb 100644
--- a/onnxcustom/training/sgd_learning_rate.py
+++ b/onnxcustom/training/sgd_learning_rate.py
@@ -10,7 +10,7 @@
     OrtValue as C_OrtValue)
 from ..utils.onnx_function import function_onnx_graph
 from ..utils.onnxruntime_helper import device_to_providers
-from .base_onnx_function import BaseLearningOnnx
+from ._base_onnx_function import BaseLearningOnnx
 
 
 class BaseLearningRate(BaseLearningOnnx):
diff --git a/onnxcustom/utils/onnx_function.py b/onnxcustom/utils/onnx_function.py
index 7306143f..80f0a81d 100644
--- a/onnxcustom/utils/onnx_function.py
+++ b/onnxcustom/utils/onnx_function.py
@@ -714,3 +714,88 @@ def _onnx_update_penalty_elastic_error(target_opset=None, dtype=numpy.float32,
         varsx, outputs=[('Y', var_type())],
         target_opset=target_opset)
     return onx
+
+
+def _onnx_grad_sigmoid_neg_log_loss_error(target_opset=None,
+                                          dtype=numpy.float32,
+                                          eps=1e-5,
+                                          weight_name=None):
+    """
+    The function the raw scores from a classifier, uses the
+    sigmoid function to compute probabilities, then the log function
+    to compute the loss. It creates the ONNX graph for this function
+    and the associated gradient of the loss against the raw scores.
+
+    Probabilites (class 1): :math:`p(s) = \\frac{1}{1 + \\exp(-s)}`.
+    Loss (for two classes): :math:`L(y, s) = (1 - y)\\log(1 - p(s)) +
+    y \\log(p(s))`.
+    Gradient :math:`\\frac{dL(y, s)}{ds} = y - p(s)`.
+    To avoid nan values, probabilies are clipped:
+    :math:`p(s) = \\max(\\min(p(s), 1 - \\epsilon), \\epsilon)`.
+    :math:`y \\in \\{0, 1\\}` (integer). *s* is a float.
+
+    :param eps: to clip probabilities and avoid computing `log(0)`
+
+    .. gdot::
+        :script: DOT-SECTION
+
+        from mlprodict.onnxrt import OnnxInference
+        from onnxcustom.utils.onnx_function import function_onnx_graph
+
+        model_onnx = function_onnx_graph('grad_sigmoid_log_loss_error')
+        oinf = OnnxInference(model_onnx, inplace=False)
+
+        print("DOT-SECTION", oinf.to_dot())
+    """
+    from onnx.mapping import NP_TYPE_TO_TENSOR_TYPE
+    from skl2onnx.algebra.onnx_ops import (
+        OnnxSub, OnnxMul, OnnxSigmoid, OnnxLog, OnnxNeg,
+        OnnxReduceMean, OnnxReshape, OnnxAdd, OnnxCast, OnnxClip)
+
+    p1c = OnnxSigmoid('X2', op_version=target_opset)
+    p1 = OnnxClip(p1c, numpy.array([eps], dtype=dtype),
+                  numpy.array([1 - eps], dtype=dtype),
+                  op_version=target_opset)
+    p0 = OnnxSub(numpy.array([1], dtype=dtype), p1,
+                 op_version=target_opset)
+    y1 = OnnxCast('X1', to=NP_TYPE_TO_TENSOR_TYPE[numpy.dtype(dtype)],
+                  op_version=target_opset)
+    y0 = OnnxSub(numpy.array([1], dtype=dtype), y1,
+                 op_version=target_opset)
+    loss_obs = OnnxAdd(
+        OnnxMul(y0, OnnxLog(p0, op_version=target_opset),
+                op_version=target_opset),
+        OnnxMul(y1, OnnxLog(p1, op_version=target_opset),
+                op_version=target_opset),
+        op_version=target_opset)
+
+    if weight_name is None:
+        loss = OnnxReduceMean(loss_obs, op_version=target_opset)
+        grad = OnnxSub(y1, p1, op_version=target_opset,
+                       output_names=['Z'])
+    else:
+        loss = OnnxReduceMean(
+            OnnxMul(loss_obs, weight_name, op_version=target_opset),
+            op_version=target_opset)
+        grad = OnnxMul(
+            OnnxSub(y1, p1, op_version=target_opset),
+            weight_name, output_names=['Z'], op_version=target_opset)
+
+    loss_neg = OnnxNeg(loss, op_version=target_opset)
+    res = OnnxReshape(loss_neg, numpy.array([-1], numpy.int64),
+                      op_version=target_opset,
+                      output_names=['Y'])
+
+    var_type_int64 = dtype_to_var_type(numpy.int64)
+    var_type = dtype_to_var_type(dtype)
+    varsx = [('X1', var_type_int64([None, None])),
+             ('X2', var_type([None, None]))]
+    if weight_name is not None:
+        varsx.append((weight_name, var_type([None])))
+    onx = res.to_onnx(
+        varsx, outputs=[('Y', var_type()), ('Z', var_type())],
+        target_opset=target_opset, other_outputs=[grad])
+    if weight_name is not None:
+        onx = add_initializer(
+            onx, weight_name, numpy.array([1], dtype=dtype))
+    return onx
diff --git a/onnxcustom/utils/onnx_helper.py b/onnxcustom/utils/onnx_helper.py
index 872c1410..1959fa50 100644
--- a/onnxcustom/utils/onnx_helper.py
+++ b/onnxcustom/utils/onnx_helper.py
@@ -10,7 +10,7 @@
 
 def onnx_rename_weights(onx):
     """
-    Renames ONNX initialiers to make sure their name
+    Renames ONNX initializers to make sure their name
     follows the alphabetical order. The model is
     modified inplace. This function calls
     :func:`onnx_rename_names
@@ -73,11 +73,17 @@ def dtype_to_var_type(dtype):
     """
     Converts a numpy dtype into a var type.
     """
-    from skl2onnx.common.data_types import FloatTensorType, DoubleTensorType
+    from skl2onnx.common.data_types import (
+        FloatTensorType, DoubleTensorType,
+        Int32TensorType, Int64TensorType)
     if dtype == numpy.float32:
         return FloatTensorType
     if dtype == numpy.float64:
         return DoubleTensorType
+    if dtype == numpy.int64:
+        return Int64TensorType
+    if dtype == numpy.int32:
+        return Int32TensorType
     raise ValueError(
         "Unexpected value dtype=%r." % dtype)
 
diff --git a/onnxcustom/utils/onnx_rewriter.py b/onnxcustom/utils/onnx_rewriter.py
new file mode 100644
index 00000000..f80568ca
--- /dev/null
+++ b/onnxcustom/utils/onnx_rewriter.py
@@ -0,0 +1,166 @@
+"""
+@file
+@brief Rewrites operator in ONNX graph.
+"""
+from onnx.helper import make_graph
+from onnx import NodeProto
+from onnx.numpy_helper import to_array, from_array
+
+
+def _unique_name(existing_names, name):
+    """
+    Returns a name different from any name in *existing_names*.
+
+    :param existing_names: set of names
+    :param name: current
+    :return: unique name
+    """
+    if name not in existing_names:
+        existing_names.add(name)
+        return name
+    name0 = name
+    i = 2
+    while name in existing_names:
+        name = "%s_%d" % (name0, i)
+        i += 1
+    existing_names.add(name)
+    return name
+
+
+def _onnx_rewrite_operator_node(existing_names, node, sub_onx):
+    """
+    Replaces a node by a subgraph.
+
+    :param existing_names: existing results names
+    :param node: onnx node to replace
+    :param sub_onx: onnx sub_graph to use as a replacement
+    :return: new_initializer, new_nodes
+    """
+    if len(node.input) != len(sub_onx.graph.input):
+        raise ValueError(
+            "Mismatch with the number of inputs for operator type %r. "
+            "%d != %d." % (
+                node.op_type, len(node.input), len(sub_onx.graph.nput)))
+    if len(node.output) != len(sub_onx.graph.output):
+        raise ValueError(
+            "Mismatch with the number of outputs for operator type %r. "
+            "%d != %d." % (
+                node.op_type, len(node.output), len(sub_onx.graph.output)))
+    replaces = {}
+    for inp, name in zip(sub_onx.graph.input, node.input):
+        replaces[inp.name] = name
+    for inp, name in zip(sub_onx.graph.output, node.output):
+        replaces[inp.name] = name
+
+    new_inits = []
+    for init in sub_onx.graph.initializer:
+        name = _unique_name(existing_names, init.name)
+        replaces[init.name] = name
+        tensor = from_array(to_array(init), name=name)
+        new_inits.append(tensor)
+
+    new_nodes = []
+    for n in sub_onx.graph.node:
+        new_node = NodeProto()
+        new_node.op_type = n.op_type
+        new_node.attribute.extend(n.attribute)  # pylint: disable=E1101
+        new_node.input.extend(  # pylint: disable=E1101
+            [replaces[i] for i in n.input])  # pylint: disable=E1101
+        new_node.domain = n.domain
+        new_out = []
+        for o in n.output:
+            if o in replaces:
+                new_out.append(replaces[o])
+            else:
+                n = _unique_name(existing_names, o)
+                new_out.append(n)
+        new_node.output.extend(new_out)  # pylint: disable=E1101
+        new_nodes.append(new_node)
+
+    return new_inits, new_nodes
+
+
+def onnx_rewrite_operator(onx, op_type, sub_onx, recursive=True, debug_info=None):
+    """
+    Replaces one operator by an onnx graph.
+
+    :param onx: onnx graph
+    :param op_type: operator type
+    :param sub_onx: onnx graph
+    :param recursive: looks into subgraphs
+    :param debug_info: unused
+    :return: modified onnx graph
+
+    .. runpython::
+        :showcode:
+
+        import numpy
+        from skl2onnx.common.data_types import FloatTensorType
+        from skl2onnx.algebra.onnx_ops import (  # pylint: disable=E0611
+            OnnxReciprocal, OnnxDiv)
+        from mlprodict.plotting.text_plot import onnx_simple_text_plot
+        from onnxcustom import get_max_opset
+        from onnxcustom.utils.onnx_rewriter import onnx_rewrite_operator
+
+        # first graph: it contains the node to replace
+        opset = get_max_opset()
+        node1 = OnnxReciprocal('X', output_names=['Y'],
+                               op_version=opset)
+        onx1 = node1.to_onnx(
+            inputs={'X': FloatTensorType()},
+            outputs={'Y': FloatTensorType()},
+            target_opset=opset)
+
+        # second graph: it contains the replacement graph
+        node2 = OnnxDiv(numpy.array([1], dtype=numpy.float32),
+                        'X', output_names=['Y'],
+                        op_version=opset)
+        onx2 = node2.to_onnx(
+            inputs={'X': FloatTensorType()},
+            outputs={'Y': FloatTensorType()},
+            target_opset=opset)
+
+        # third graph: the modified graph
+        onx3 = onnx_rewrite_operator(onx1, 'Reciprocal', onx2)
+        print(onnx_simple_text_plot(onx3))
+    """
+    from mlprodict.onnx_tools.optim._onnx_optimisation_common import (  # pylint: disable=C0415
+        _apply_remove_node_fct_node, _apply_optimisation_on_graph)
+
+    if hasattr(onx, 'graph'):
+        fct = (lambda graph, recursive=False, debug_info=None:
+               onnx_rewrite_operator(
+                   graph, op_type, sub_onx, recursive=recursive,
+                   debug_info=debug_info))
+        return _apply_optimisation_on_graph(fct, onx, recursive=recursive)
+
+    existing_names = set()
+    for node in onx.node:
+        existing_names.update(node.input)
+        existing_names.update(node.output)
+
+    nodes = list(onx.node)
+    new_nodes = []
+    new_inits = list(onx.initializer)
+    for i, node in enumerate(nodes):
+        if node.op_type != op_type:
+            new_nodes.append(node)
+            continue
+        inits, newn = _onnx_rewrite_operator_node(
+            existing_names, node, sub_onx)
+        new_inits.extend(inits)
+        new_nodes.extend(newn)
+
+    if recursive:
+        # Handles subgraphs.
+        for i in range(len(new_nodes)):  # pylint: disable=C0200
+            node = nodes[i]
+            if node is None or not (node.attribute):  # pylint: disable=C0325
+                continue
+            nodes[i] = _apply_remove_node_fct_node(
+                onnx_rewrite_operator, node, recursive=True,
+                debug_info=None)
+
+    graph = make_graph(
+        new_nodes, onx.name, onx.input, onx.output, new_inits)
+    return graph
diff --git a/onnxcustom/utils/orttraining_helper.py b/onnxcustom/utils/orttraining_helper.py
index 74441d76..c4f9c437 100644
--- a/onnxcustom/utils/orttraining_helper.py
+++ b/onnxcustom/utils/orttraining_helper.py
@@ -13,13 +13,12 @@
 from onnx import TensorProto
 
 
-def _unique_name(existing_names, name, add=True):
+def _unique_name(existing_names, name):
     """
     Returns a name different from any name in *existing_names*.
 
     :param existing_names: set of names
     :param name: current
-    :param add: add the name of the list of existing names
     :return: unique name
     """
     if name not in existing_names:
@@ -134,6 +133,85 @@ def _loss_elastic(existing_names, elem, shape,
         [make_tensor_value_info(loss_name, elem, [1, 1])])
 
 
+def _loss_log(existing_names, elem, shape,
+              output_name, label_name,
+              weight_name, loss_name,
+              eps=1e-6):
+    """
+    This only works for a binary classification.
+    The log loss is `'log(yt, yp) = (1-yt)\\log(1-yp) - yt\\log(yp)`,
+    this only works for a binary classification where *yp* is the
+    predicted probability, *yt* is the expected probability.
+    *yt* is expected to be binary, *yp* is a matrix with two
+    columns, the sum on every line is 1.
+    Parameter *eps* is used to avoid computing *log(0)*.
+    """
+    if output_name == 'output_label':
+        raise RuntimeError(
+            "output_name=%r, log loss does not work on labels."
+            "" % output_name)
+    dtype = TENSOR_TYPE_TO_NP_TYPE[elem]
+    one_name = _unique_name(existing_names, "one_name")
+    eps_name = _unique_name(existing_names, "eps_name")
+    eps1_name = _unique_name(existing_names, "eps1_name")
+    axes_name = _unique_name(existing_names, "axes_name")
+
+    eps_init = from_array(numpy.array([eps], dtype=dtype), name=eps_name)
+    one_init = from_array(numpy.array([1], dtype=dtype), name=one_name)
+    eps1_init = from_array(
+        numpy.array([1 - eps], dtype=dtype), name=eps1_name)
+    axes_init = from_array(
+        numpy.array([1], dtype=numpy.int64), name=axes_name)
+
+    clip_name = _unique_name(existing_names, "clip_name")
+    clip_red_name = _unique_name(existing_names, "clip_red_name")
+    new_output_name = _unique_name(existing_names, "new_output_name")
+    cast_name = _unique_name(existing_names, "cast_name")
+    log_name = _unique_name(existing_names, "log_name")
+    subl_name = _unique_name(existing_names, "subl_name")
+    conc_name = _unique_name(existing_names, "conc_name")
+    mul_name = _unique_name(existing_names, "mul_name")
+    like_name = _unique_name(existing_names, "like_name")
+
+    nodes = [
+        make_node(
+            'Clip', [output_name, eps_name, eps1_name], [clip_name]),
+        make_node(
+            'ReduceSum', [clip_name, axes_name], [clip_red_name], keepdims=1),
+        make_node('Div', [clip_name, clip_red_name], [new_output_name]),
+        make_node('Log', [new_output_name], [log_name]),
+        make_node('Cast', [label_name], [cast_name], to=elem),
+        make_node('Sub', [one_name, cast_name], [subl_name]),
+        make_node('Concat', [subl_name, cast_name], [conc_name], axis=1),
+        make_node('Mul', [log_name, conc_name], [mul_name]),
+        make_node(
+            'ReduceSum', [mul_name, axes_name], [like_name], keepdims=1)]
+
+    inputs = [make_tensor_value_info(label_name, TensorProto.INT64, shape)]
+
+    if weight_name is not None:
+        inputs.append(
+            make_tensor_value_info(weight_name, elem, [shape[0]]))
+        likew_name = _unique_name(existing_names, "likew_name")
+        nodes.append(
+            make_node('Mul', [like_name, weight_name], [likew_name]))
+        like_name = likew_name
+
+    shape_name = _unique_name(existing_names, "shape_name")
+    onx_shape = from_array(
+        numpy.array([1, 1], dtype=numpy.int64), name=shape_name)
+    reduced_loss = _unique_name(existing_names, "reduced_loss")
+    neg_reduced_loss = _unique_name(existing_names, "neg_reduced_loss")
+    nodes.extend([
+        make_node('ReduceMean', [like_name], [reduced_loss]),
+        make_node('Neg', [reduced_loss], [neg_reduced_loss]),
+        make_node('Reshape', [neg_reduced_loss, shape_name], [loss_name])])
+
+    return (
+        [eps_init, eps1_init, one_init, axes_init, onx_shape],
+        inputs, nodes, [make_tensor_value_info(loss_name, elem, [1, 1])])
+
+
 def penalty_loss_onnx(name, dtype, l1=None, l2=None, existing_names=None):
     """
     Returns onnx nodes to compute
@@ -210,10 +288,59 @@ def penalty_loss_onnx(name, dtype, l1=None, l2=None, existing_names=None):
     return inits, nodes
 
 
+def get_train_initializer(onx):
+    """
+    Returns the list of initializers to train.
+
+    :return: dictionary `{name: (value, tensor)}`
+
+    The function walk through the list of initializers and
+    returns all tensors with elements from types float or double.
+    """
+    res = OrderedDict()
+    for init in onx.graph.initializer:
+        if init.data_type in (
+                TensorProto.FLOAT16,  # pylint: disable=E1101
+                TensorProto.FLOAT,  # pylint: disable=E1101
+                TensorProto.DOUBLE):  # pylint: disable=E1101
+            res[init.name] = (to_array(init), init)
+    return res
+
+
+def _rewrite_op_no_grad(onx):
+    """
+    Rewrites operators with no gradient.
+    """
+    set_types = set(n.name for n in onx.graph.node)
+    if "Reciprocal" in set_types:
+        from skl2onnx.algebra.onnx_ops import OnnxDiv  # pylint: disable=E0611
+        from skl2onnx.common.data_types import FloatTensorType
+        from .onnx_rewriter import onnx_rewrite_operator
+
+        opset = None
+        for op in onx.opset_import:
+            if op.domain in ('', 'ai.onnx'):
+                opset = op.version
+        if opset is None:
+            from .. import get_max_opset
+            opset = get_max_opset()
+
+        node = OnnxDiv(numpy.array([1], dtype=numpy.float32),
+                       'X', output_names=['Y'],
+                       op_version=opset)
+        rewrite_onx = node.to_onnx(
+            inputs={'X': FloatTensorType()},
+            outputs={'Y': FloatTensorType()},
+            target_opset=opset)
+        onx = onnx_rewrite_operator(onx, 'Reciprocal', rewrite_onx)
+
+    return onx
+
+
 def add_loss_output(onx, score_name='squared_error',
                     loss_name='loss', label_name='label',
-                    weight_name=None,
-                    penalty=None, **kwargs):
+                    weight_name=None, penalty=None,
+                    output_index=None, **kwargs):
     """
     Modifies an ONNX graph to add operators to score and allow training.
 
@@ -228,6 +355,10 @@ def add_loss_output(onx, score_name='squared_error',
         or `{ weight_name: beta}`,
         it adds a L1 and/or L2 penalty to one input or initializer,
         penalty = :math:`|w| \\alpha + w^2 \\beta`
+    :param output_index: the output used to compute the loss,
+        if None, the function assumes there is only one output,
+        it must be specified if there are more than 1,
+        it can be an integer or a string (output name)
     :param kwargs: additional arguments for losses (see below)
     :return: modified graph
 
@@ -241,6 +372,11 @@ def add_loss_output(onx, score_name='squared_error',
       is not None
     * `'elastic'`: mixture of losses, kwargs must define
       *l1_weight* and *l2_weight*, undefined, default value are 0.5
+    * `'log(yt, yp)'`: log loss :math:`(1-yt)\\log(1-yp) - yt\\log(yp)`,
+        this only works for a binary classification where *yp* is the
+        predicted probability, *yt* is the expected probability.
+        *yt* is expected to be binary, *yp* is a matrix with two
+        columns, the sum on every line is 1.
 
     See example :ref:`l-orttraining-nn-gpu`.
     Next example shows the loss with L1 and L2 loss.
@@ -311,11 +447,45 @@ def add_loss_output(onx, score_name='squared_error',
 
         print("DOT-SECTION", OnnxInference(onx_loss).to_dot())
     """
-    outputs = onx.graph.output
-    if len(outputs) != 1:
-        raise ValueError(  # pragma: no cover
-            "Unable to guess the output to compare to the "
-            "expacted labels among %r." % (o.name for o in outputs))
+    from mlprodict.onnx_tools.optim import onnx_remove_node_unused
+
+    # rename every intermediate output call label
+    def _replace(ens):
+        for i in range(len(ens)):  # pylint: disable=C0200
+            if ens[i] == 'label':
+                ens[i] = '_label_'
+
+    for node in onx.graph.node:
+        if "_label_" in node.input or "_label_" in node.output:
+            raise RuntimeError(  # pragma: no cover
+                "One intermediate result contains '_label_'. "
+                "It should be removed manually.\n%r" % node)
+        _replace(node.input)
+        _replace(node.output)
+
+    if output_index is None:
+        if len(onx.graph.output) != 1:
+            raise ValueError(  # pragma: no cover
+                "Unable to guess the output to compare to the "
+                "expacted labels among %r." % (
+                    [o.name for o in onx.graph.output]))
+        outputs = onx.graph.output
+        output_index = 0
+    elif isinstance(output_index, int):
+        outputs = [onx.graph.output[output_index]]
+    elif isinstance(output_index, str):
+        outputs = [(i, o) for i, o in enumerate(onx.graph.output)
+                   if o.name == output_index]
+        if len(outputs) != 1:
+            raise ValueError(  # pragma: no cover
+                "Unable to find output %r in %r." % (
+                    output_index, [o.name for o in onx.graph.output]))
+        output_index = outputs[0][0]
+        outputs = [outputs[0][1]]
+    else:
+        raise TypeError(
+            "output_index must be an integer or a str not %r."
+            "" % type(output_index))
 
     existing_names = []
     for node in onx.graph.node:
@@ -323,10 +493,15 @@ def add_loss_output(onx, score_name='squared_error',
         existing_names.extend(node.input)
     existing_names = set(existing_names)
 
-    output_name = onx.graph.output[0].name
-    elem = onx.graph.output[0].type.tensor_type.elem_type
+    output_onx = onx.graph.output[output_index]
+    output_name = output_onx.name
+    elem = output_onx.type.tensor_type.elem_type
+    if elem == 0:
+        raise TypeError(  # pragma: no cover
+            "Unable to guess inut tensor type from %r."
+            "" % output_onx)
     shape = []
-    for d in onx.graph.output[0].type.tensor_type.shape.dim:
+    for d in output_onx.type.tensor_type.shape.dim:
         shape.append(d.dim_value if d.dim_value > 0 else None)
 
     if score_name in ('squared_error', 'l2'):
@@ -341,6 +516,11 @@ def add_loss_output(onx, score_name='squared_error',
         inits, inputs, nodes, outputs = _loss_elastic(
             existing_names, elem, shape, output_name, label_name,
             weight_name, loss_name, **kwargs)
+    elif score_name == 'log':
+        shape = (None, 1)
+        inits, inputs, nodes, outputs = _loss_log(
+            existing_names, elem, shape, output_name, label_name,
+            weight_name, loss_name, **kwargs)
     else:
         raise NotImplementedError(  # pragma: no cover
             "Unexpected %r value for score_name." % score_name)
@@ -388,7 +568,7 @@ def add_loss_output(onx, score_name='squared_error',
         list(onx.graph.node) + nodes,
         onx.graph.name,
         list(onx.graph.input) + inputs,
-        outputs + list(onx.graph.output),
+        outputs + [onx.graph.output[output_index]],
         inits)
     onnx_model = make_model(graph)
     onnx_model.ir_version = onx.ir_version
@@ -407,23 +587,4 @@ def add_loss_output(onx, score_name='squared_error',
         op_set = onnx_model.opset_import.add()  # pylint: disable=E1101
         op_set.domain = oimp.domain
         op_set.version = oimp.version
-    return onnx_model
-
-
-def get_train_initializer(onx):
-    """
-    Returns the list of initializers to train.
-
-    :return: dictionary `{name: (value, tensor)}`
-
-    The function walk through the list of initializers and
-    returns all tensors with elements from types float or double.
-    """
-    res = OrderedDict()
-    for init in onx.graph.initializer:
-        if init.data_type in (
-                TensorProto.FLOAT16,  # pylint: disable=E1101
-                TensorProto.FLOAT,  # pylint: disable=E1101
-                TensorProto.DOUBLE):  # pylint: disable=E1101
-            res[init.name] = (to_array(init), init)
-    return res
+    return _rewrite_op_no_grad(onnx_remove_node_unused(onnx_model))


Test patch:
diff --git a/_unittests/ut_module/test_onnx_runtimes.py b/_unittests/ut_module/test_onnx_runtimes.py
new file mode 100644
index 00000000..a8336ccc
--- /dev/null
+++ b/_unittests/ut_module/test_onnx_runtimes.py
@@ -0,0 +1,42 @@
+"""
+@brief      test log(time=0s)
+"""
+import unittest
+import numpy
+from scipy.special import expit  # pylint: disable=E0611
+from pyquickhelper.pycode import ExtTestCase
+from skl2onnx.algebra.onnx_ops import OnnxSigmoid  # pylint: disable=E0611
+from skl2onnx.common.data_types import FloatTensorType
+from mlprodict.onnxrt import OnnxInference
+from onnxcustom import get_max_opset
+
+
+class TestOnnxRuntimes(ExtTestCase):
+    """Test style."""
+
+    def test_check(self):
+        opset = get_max_opset()
+        min_values = [-41.621277, -40.621277, -30.621277, -20.621277,
+                      -19, -18, -17, -15, -14, -13, -12, -11, -10, -5, -2]
+        data = numpy.array(
+            [[0]],
+            dtype=numpy.float32)
+
+        node = OnnxSigmoid('X', op_version=opset, output_names=['Y'])
+        onx = node.to_onnx({'X': FloatTensorType()},
+                           {'Y': FloatTensorType()},
+                           target_opset=opset)
+        rts = ['numpy', 'python', 'onnxruntime1']
+        for mv in min_values:
+            data[:, 0] = mv
+            for rt in rts:
+                if rt == 'numpy':
+                    y = expit(data)
+                else:
+                    oinf = OnnxInference(onx, runtime=rt)
+                    y = oinf.run({'X': data})['Y']
+                self.assertNotEmpty(y)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/_unittests/ut_training/test_optimizers.py b/_unittests/ut_training/test_optimizers.py
index d441aaa6..140573b5 100644
--- a/_unittests/ut_training/test_optimizers.py
+++ b/_unittests/ut_training/test_optimizers.py
@@ -1,5 +1,5 @@
 """
-@brief      test log(time=8s)
+@brief      test log(time=9s)
 """
 import os
 import unittest
diff --git a/_unittests/ut_training/test_optimizers_classification.py b/_unittests/ut_training/test_optimizers_classification.py
new file mode 100644
index 00000000..0bf9cfc0
--- /dev/null
+++ b/_unittests/ut_training/test_optimizers_classification.py
@@ -0,0 +1,95 @@
+"""
+@brief      test log(time=8s)
+"""
+import unittest
+from pyquickhelper.pycode import ExtTestCase, get_temp_folder
+import numpy
+from onnx.helper import set_model_props
+from sklearn.datasets import make_classification
+from sklearn.model_selection import train_test_split
+from sklearn.linear_model import SGDClassifier
+from mlprodict.onnx_conv import to_onnx
+from mlprodict.plotting.text_plot import onnx_simple_text_plot
+from mlprodict.onnx_tools.onnx_manipulations import select_model_inputs_outputs
+# from mlprodict.onnxrt import OnnxInference
+from onnxcustom import __max_supported_opset__ as opset
+try:
+    from onnxruntime import TrainingSession
+except ImportError:
+    # onnxruntime not training
+    TrainingSession = None
+
+
+class TestOptimizersClassification(ExtTestCase):
+
+    @unittest.skipIf(TrainingSession is None, reason="not training")
+    def test_ort_gradient_optimizers_binary(self):
+        from onnxcustom.utils.orttraining_helper import add_loss_output
+        from onnxcustom.training.optimizers import OrtGradientOptimizer
+        X, y = make_classification(  # pylint: disable=W0632
+            100, n_features=10, random_state=0)
+        X = X.astype(numpy.float32)
+        y = y.astype(numpy.int64)
+        X_train, _, y_train, __ = train_test_split(X, y)
+        reg = SGDClassifier(loss='log')
+        reg.fit(X_train, y_train)
+        onx = to_onnx(reg, X_train, target_opset=opset,
+                      black_op={'LinearClassifier'},
+                      options={'zipmap': False})
+        set_model_props(onx, {'info': 'unit test'})
+        onx_loss = add_loss_output(onx, 'log', output_index=1)
+        inits = ['intercept', 'coef']
+        train_session = OrtGradientOptimizer(
+            onx_loss, inits, learning_rate=1e-3)
+        self.assertRaise(lambda: train_session.get_state(), AttributeError)
+        train_session.fit(X_train, y_train.reshape((-1, 1)), use_numpy=True)
+        state_tensors = train_session.get_state()
+        self.assertEqual(len(state_tensors), 2)
+        r = repr(train_session)
+        self.assertIn("OrtGradientOptimizer(model_onnx=", r)
+        self.assertIn("learning_rate='invscaling'", r)
+        losses = train_session.train_losses_
+        self.assertGreater(len(losses), 1)
+        self.assertFalse(any(map(numpy.isnan, losses)))
+
+    @unittest.skipIf(TrainingSession is None, reason="not training")
+    def test_ort_gradient_optimizers_fw_nesterov_binary(self):
+        from onnxcustom.training.optimizers_partial import (
+            OrtGradientForwardBackwardOptimizer)
+        from onnxcustom.training.sgd_learning_rate import (
+            LearningRateSGDNesterov)
+        from onnxcustom.training.sgd_learning_loss import NegLogLearningLoss
+        X, y = make_classification(  # pylint: disable=W0632
+            100, n_features=10, random_state=0)
+        X = X.astype(numpy.float32)
+        y = y.astype(numpy.int64)
+        X_train, _, y_train, __ = train_test_split(X, y)
+        reg = SGDClassifier(loss='log')
+        reg.fit(X_train, y_train)
+        onx = to_onnx(reg, X_train, target_opset=opset,
+                      black_op={'LinearRegressor'},
+                      options={'zipmap': False,
+                               'raw_scores': True})
+        onx = select_model_inputs_outputs(onx, outputs=['score'])
+        self.assertIn("output: name='score'",
+                      onnx_simple_text_plot(onx))
+        set_model_props(onx, {'info': 'unit test'})
+        inits = ['coef', 'intercept']
+
+        train_session = OrtGradientForwardBackwardOptimizer(
+            onx, inits,
+            learning_rate=LearningRateSGDNesterov(
+                1e-4, nesterov=False, momentum=0.9),
+            learning_loss=NegLogLearningLoss(),
+            warm_start=False, max_iter=100, batch_size=10)
+        self.assertIsInstance(train_session.learning_loss, NegLogLearningLoss)
+        self.assertEqual(train_session.learning_loss.eps, 1e-5)
+        train_session.fit(X, y)
+        temp = get_temp_folder(
+            __file__, "temp_ort_gradient_optimizers_fw_nesterov_binary")
+        train_session.save_onnx_graph(temp)
+
+
+if __name__ == "__main__":
+    TestOptimizersClassification().test_ort_gradient_optimizers_fw_nesterov_binary()
+    unittest.main()
diff --git a/_unittests/ut_training/test_orttraining_forward_backward.py b/_unittests/ut_training/test_orttraining_forward_backward.py
index 4ecff37c..cbc57faf 100644
--- a/_unittests/ut_training/test_orttraining_forward_backward.py
+++ b/_unittests/ut_training/test_orttraining_forward_backward.py
@@ -401,4 +401,4 @@ def test_forward_training_logreg(self):
 
 if __name__ == "__main__":
     # TestOrtTrainingForwardBackward().forward_no_training(verbose=True)
-    unittest.main(verbosity=2)
+    unittest.main(verbosity=2, failfast=True)
diff --git a/_unittests/ut_utils/test_onnx_function.py b/_unittests/ut_utils/test_onnx_function.py
index 995cc9aa..62b4626e 100644
--- a/_unittests/ut_utils/test_onnx_function.py
+++ b/_unittests/ut_utils/test_onnx_function.py
@@ -3,6 +3,7 @@
 """
 import unittest
 import numpy
+from scipy.special import expit  # pylint: disable=E0611
 from onnxruntime import InferenceSession, SessionOptions
 from pyquickhelper.pycode import ExtTestCase
 from mlprodict.onnxrt import OnnxInference
@@ -108,13 +109,30 @@ def common_check_alpha(self, name, fct):
     def test_grad_onnx_axpy(self):
         self.common_check_alpha("axpy", lambda x1, x2, alpha: x1 * alpha + x2)
 
-    def common_check_2(self, name, fct, weight_name=None, **kwargs):
+    def common_check_2(self, name, fct, weight_name=None,
+                       verbose=0, classification=False, **kwargs):
         onx = function_onnx_graph(
             name, target_opset=get_max_opset(),
             dtype=numpy.float32, weight_name=weight_name,
             **kwargs)
-        x1 = numpy.random.randn(10, 1).astype(numpy.float32)
-        x2 = numpy.random.randn(10, 1).astype(numpy.float32)
+        if verbose > 0:
+            with open(name + ".onnx", "wb") as f:
+                f.write(onx.SerializeToString())
+        if classification:
+            N = 10
+            p = numpy.random.randn(N, 1).astype(numpy.float32)
+            p[0, :] = 0
+            p[1, :] = 100
+            p[2, :] = -100
+            p[3, :] = 1
+            p[4, :] = -1
+            y = (numpy.random.randn(N, 1).astype(numpy.float32) > 0).astype(
+                numpy.int64)
+            x2 = p
+            x1 = y
+        else:
+            x1 = numpy.random.randn(10, 1).astype(numpy.float32)
+            x2 = numpy.random.randn(10, 1).astype(numpy.float32)
         w = numpy.random.rand(10).astype(numpy.float32)
         if weight_name is None:
             exp_loss, exp_grad = fct(x1, x2)
@@ -122,12 +140,13 @@ def common_check_2(self, name, fct, weight_name=None, **kwargs):
             exp_loss, exp_grad = fct(x1, x2, w)
 
         oinf = OnnxInference(onx)
+        run_params = dict(verbose=verbose, fLOG=print) if verbose > 0 else {}
         if weight_name is None:
-            got = oinf.run({'X1': x1, 'X2': x2})
+            got = oinf.run({'X1': x1, 'X2': x2}, **run_params)
         else:
-            got = oinf.run({'X1': x1, 'X2': x2, 'weight': w})
-        self.assertEqualArray(exp_loss, got['Y'], decimal=5)
+            got = oinf.run({'X1': x1, 'X2': x2, 'weight': w}, **run_params)
         self.assertEqualArray(exp_grad, got['Z'], decimal=5)
+        self.assertEqualArray(exp_loss, got['Y'], decimal=5)
 
         providers = device_to_providers('cpu')
         so = SessionOptions()
@@ -384,6 +403,19 @@ def fct(x):
         got = sess.run(None, {'X': x})
         self.assertEqualArray(exp_loss, got[0], decimal=5)
 
+    def test_grad_sigmoid_neg_log_loss_error(self):
+
+        def loss(x1, x2, eps=1e-5):
+            pr = expit(x2)
+            cl = numpy.clip(pr, eps, 1 - eps)
+            lo = - (1 - x1) * numpy.log(1 - cl) - x1 * numpy.log(cl)
+            return lo
+
+        self.common_check_2(
+            "grad_sigmoid_neg_log_loss_error",
+            lambda x1, x2: (loss(x1, x2).mean(), x1 - expit(x2)),
+            classification=True)
+
 
 if __name__ == "__main__":
     unittest.main()
diff --git a/_unittests/ut_utils/test_onnx_helper.py b/_unittests/ut_utils/test_onnx_helper.py
index e16774a8..b1c0bee5 100644
--- a/_unittests/ut_utils/test_onnx_helper.py
+++ b/_unittests/ut_utils/test_onnx_helper.py
@@ -4,7 +4,8 @@
 import unittest
 import numpy
 from pyquickhelper.pycode import ExtTestCase
-from skl2onnx.common.data_types import FloatTensorType, DoubleTensorType
+from skl2onnx.common.data_types import (
+    FloatTensorType, DoubleTensorType, Int64TensorType)
 from skl2onnx.algebra.onnx_ops import (  # pylint: disable=E0611
     OnnxRelu, OnnxMatMul)
 from onnxcustom.utils.onnx_helper import (
@@ -36,8 +37,9 @@ def test_onnx_rename_weights(self):
     def test_dtype_to_var_type(self):
         self.assertEqual(dtype_to_var_type(numpy.float32), FloatTensorType)
         self.assertEqual(dtype_to_var_type(numpy.float64), DoubleTensorType)
-        self.assertRaise(lambda: dtype_to_var_type(numpy.int64), ValueError)
+        self.assertEqual(dtype_to_var_type(numpy.int64), Int64TensorType)
         self.assertEqual(proto_type_to_dtype('tensor(double)'), numpy.float64)
+        self.assertRaise(lambda: dtype_to_var_type(numpy.int8), ValueError)
 
     def test_proto_type_to_dtype(self):
         self.assertEqual(proto_type_to_dtype(1), numpy.float32)
diff --git a/_unittests/ut_utils/test_onnx_orttraining.py b/_unittests/ut_utils/test_onnx_orttraining.py
deleted file mode 100644
index 7f73a486..00000000
--- a/_unittests/ut_utils/test_onnx_orttraining.py
+++ /dev/null
@@ -1,62 +0,0 @@
-"""
-@brief      test log(time=9s)
-"""
-
-import unittest
-from pyquickhelper.pycode import ExtTestCase
-import numpy
-from sklearn.datasets import make_regression
-from sklearn.model_selection import train_test_split
-from sklearn.linear_model import LinearRegression
-from sklearn.metrics import mean_squared_error
-from mlprodict.onnx_conv import to_onnx
-from mlprodict.onnxrt import OnnxInference
-from onnxcustom import __max_supported_opset__ as opset
-try:
-    from onnxruntime import TrainingSession
-except ImportError:
-    # onnxruntime not training
-    TrainingSession = None
-
-
-class TestOrtTraining(ExtTestCase):
-
-    @unittest.skipIf(TrainingSession is None, reason="not training")
-    def test_add_loss_output(self):
-        from onnxcustom.utils.orttraining_helper import add_loss_output
-        X, y = make_regression(  # pylint: disable=W0632
-            100, n_features=10, bias=2)
-        X = X.astype(numpy.float32)
-        y = y.astype(numpy.float32)
-        X_train, X_test, y_train, y_test = train_test_split(X, y)
-        reg = LinearRegression()
-        reg.fit(X_train, y_train)
-        reg.coef_ = reg.coef_.reshape((1, -1))
-        onx = to_onnx(reg, X_train, target_opset=opset,
-                      black_op={'LinearRegressor'})
-        onx_loss = add_loss_output(onx)
-        oinf = OnnxInference(onx_loss)
-        output = oinf.run({'X': X_test, 'label': y_test.reshape((-1, 1))})
-        loss = output['loss']
-        skl_loss = mean_squared_error(reg.predict(X_test), y_test)
-        self.assertLess(numpy.abs(skl_loss - loss[0, 0]), 1e-5)
-
-    @unittest.skipIf(TrainingSession is None, reason="not training")
-    def test_get_train_initializer(self):
-        from onnxcustom.utils.orttraining_helper import get_train_initializer
-        X, y = make_regression(  # pylint: disable=W0632
-            100, n_features=10, bias=2)
-        X = X.astype(numpy.float32)
-        y = y.astype(numpy.float32)
-        X_train, _, y_train, __ = train_test_split(X, y)
-        reg = LinearRegression()
-        reg.fit(X_train, y_train)
-        reg.coef_ = reg.coef_.reshape((1, -1))
-        onx = to_onnx(reg, X_train, target_opset=opset,
-                      black_op={'LinearRegressor'})
-        inits = get_train_initializer(onx)
-        self.assertEqual({'intercept', 'coef'}, set(inits))
-
-
-if __name__ == "__main__":
-    unittest.main()
diff --git a/_unittests/ut_utils/test_onnx_writer.py b/_unittests/ut_utils/test_onnx_writer.py
new file mode 100644
index 00000000..53972dd7
--- /dev/null
+++ b/_unittests/ut_utils/test_onnx_writer.py
@@ -0,0 +1,49 @@
+"""
+@brief      test log(time=9s)
+"""
+import unittest
+import numpy
+from pyquickhelper.pycode import ExtTestCase
+from skl2onnx.common.data_types import FloatTensorType
+from skl2onnx.algebra.onnx_ops import (  # pylint: disable=E0611
+    OnnxReciprocal, OnnxDiv)
+from mlprodict.onnxrt import OnnxInference
+from onnxcustom import get_max_opset
+from onnxcustom.utils.onnx_rewriter import onnx_rewrite_operator
+
+
+class TestOnnxWriter(ExtTestCase):
+
+    def test_onnx_rewrite_operator(self):
+        opset = get_max_opset()
+        node1 = OnnxReciprocal('X', output_names=['Y'],
+                               op_version=opset)
+        onx1 = node1.to_onnx(
+            inputs={'X': FloatTensorType()},
+            outputs={'Y': FloatTensorType()},
+            target_opset=opset)
+        onx1.graph.name = "jjj"
+        oinf1 = OnnxInference(onx1)
+
+        node2 = OnnxDiv(numpy.array([1], dtype=numpy.float32),
+                        'X', output_names=['Y'],
+                        op_version=opset)
+        onx2 = node2.to_onnx(
+            inputs={'X': FloatTensorType()},
+            outputs={'Y': FloatTensorType()},
+            target_opset=opset)
+        oinf2 = OnnxInference(onx2)
+        X = numpy.array([[5, 6]], dtype=numpy.float32)
+        y1 = oinf1.run({'X': X})['Y']
+        y2 = oinf2.run({'X': X})['Y']
+        self.assertEqualArray(y1, y2)
+
+        onx3 = onnx_rewrite_operator(onx1, 'Reciprocal', onx2)
+        self.assertNotIn('Reciprocal', str(onx3))
+        oinf3 = OnnxInference(onx3)
+        y3 = oinf3.run({'X': X})['Y']
+        self.assertEqualArray(y1, y3)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/_unittests/ut_utils/test_orttraining_helper.py b/_unittests/ut_utils/test_orttraining_helper.py
new file mode 100644
index 00000000..fb18f8ac
--- /dev/null
+++ b/_unittests/ut_utils/test_orttraining_helper.py
@@ -0,0 +1,110 @@
+"""
+@brief      test log(time=5s)
+"""
+import unittest
+from pyquickhelper.pycode import ExtTestCase
+import numpy
+from sklearn.datasets import make_regression, make_classification
+from sklearn.model_selection import train_test_split
+from sklearn.linear_model import LinearRegression, LogisticRegression
+from sklearn.metrics import mean_squared_error, log_loss
+from skl2onnx.algebra.onnx_ops import OnnxIdentity  # pylint: disable=E0611
+from skl2onnx.common.data_types import DoubleTensorType
+from mlprodict.onnx_conv import to_onnx
+from mlprodict.onnxrt import OnnxInference
+from mlprodict.plotting.text_plot import onnx_simple_text_plot
+from onnxcustom import __max_supported_opset__ as opset
+try:
+    from onnxruntime import TrainingSession
+except ImportError:
+    # onnxruntime not training
+    TrainingSession = None
+
+
+class TestOrtTraining(ExtTestCase):
+
+    @unittest.skipIf(TrainingSession is None, reason="not training")
+    def test_add_loss_output_reg(self):
+        from onnxcustom.utils.orttraining_helper import add_loss_output
+        X, y = make_regression(  # pylint: disable=W0632
+            100, n_features=10, bias=2)
+        X = X.astype(numpy.float32)
+        y = y.astype(numpy.float32)
+        X_train, X_test, y_train, y_test = train_test_split(X, y)
+        reg = LinearRegression()
+        reg.fit(X_train, y_train)
+        reg.coef_ = reg.coef_.reshape((1, -1))
+        onx = to_onnx(reg, X_train, target_opset=opset,
+                      black_op={'LinearRegressor'})
+        onx_loss = add_loss_output(onx)
+        oinf = OnnxInference(onx_loss)
+        output = oinf.run({'X': X_test, 'label': y_test.reshape((-1, 1))})
+        loss = output['loss']
+        skl_loss = mean_squared_error(reg.predict(X_test), y_test)
+        self.assertLess(numpy.abs(skl_loss - loss[0, 0]), 1e-5)
+
+    @unittest.skipIf(TrainingSession is None, reason="not training")
+    def test_get_train_initializer(self):
+        from onnxcustom.utils.orttraining_helper import get_train_initializer
+        X, y = make_regression(  # pylint: disable=W0632
+            100, n_features=10, bias=2)
+        X = X.astype(numpy.float32)
+        y = y.astype(numpy.float32)
+        X_train, _, y_train, __ = train_test_split(X, y)
+        reg = LinearRegression()
+        reg.fit(X_train, y_train)
+        reg.coef_ = reg.coef_.reshape((1, -1))
+        onx = to_onnx(reg, X_train, target_opset=opset,
+                      black_op={'LinearRegressor'})
+        inits = get_train_initializer(onx)
+        self.assertEqual({'intercept', 'coef'}, set(inits))
+
+    @unittest.skipIf(TrainingSession is None, reason="not training")
+    def test_add_log_loss(self):
+        from onnxcustom.utils.orttraining_helper import add_loss_output
+        ide = OnnxIdentity("X", op_version=opset, output_names=['Y'])
+        onx = ide.to_onnx(inputs={'X': DoubleTensorType()},
+                          outputs={'Y': DoubleTensorType()},
+                          target_opset=opset)
+        onx_loss = add_loss_output(onx, 'log', eps=1e-6)
+        x1 = numpy.array([0, 0, 0.2, 0.5, 0.8, 1, 1])
+        X = numpy.vstack([1 - x1, x1]).T.astype(numpy.float64)
+        y = numpy.array([0, 1, 0, 1, 1, 1, 0], dtype=numpy.int64)
+        oinf = OnnxInference(onx_loss)
+        output = oinf.run({'X': X, 'label': y.reshape((-1, 1))})
+        loss = output['loss']
+        skl_loss = log_loss(y, X[:, 1], eps=1e-6)
+        self.assertLess(numpy.abs(skl_loss - loss[0, 0]), 1e-5)
+
+    @unittest.skipIf(TrainingSession is None, reason="not training")
+    def test_add_loss_output_cls(self):
+        from onnxcustom.utils.orttraining_helper import add_loss_output
+        X, y = make_classification(  # pylint: disable=W0632
+            100, n_features=10)
+        X = X.astype(numpy.float32)
+        y = y.astype(numpy.int64)
+        X_train, X_test, y_train, y_test = train_test_split(X, y)
+        reg = LogisticRegression()
+        reg.fit(X_train, y_train)
+        reg.coef_ = reg.coef_.reshape((1, -1))
+        onx = to_onnx(reg, X_train, target_opset=opset,
+                      black_op={'LinearClassifier'},
+                      options={'zipmap': False})
+        onx_loss = add_loss_output(
+            onx, 'log', output_index='probabilities', eps=1 - 6)
+        try:
+            text = onnx_simple_text_plot(onx_loss)
+        except RuntimeError:
+            text = ""
+        if text:
+            self.assertIn("Clip(probabilities", text)
+
+        oinf = OnnxInference(onx_loss)
+        output = oinf.run({'X': X_test, 'label': y_test.reshape((-1, 1))})
+        loss = output['loss']
+        skl_loss = log_loss(y_test, reg.predict_proba(X_test), eps=1 - 6)
+        self.assertLess(numpy.abs(skl_loss - loss[0, 0]), 1e-5)
+
+
+if __name__ == "__main__":
+    unittest.main()